{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time, sys\n",
    "\n",
    "# Import the Qiskit \n",
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, QiskitError\n",
    "from qiskit import execute, IBMQ, BasicAer, transpiler, Aer\n",
    "from qiskit.providers.ibmq import least_busy\n",
    "from qiskit.tools.monitor import job_monitor\n",
    "from qiskit.mapper import Layout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define your backend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Remote backend \"ibmqx_qasm_simulator\" could not be instantiated due to an invalid config: {'conditional': ['Missing data for required field.'], 'basis_gates': ['Missing data for required field.'], 'local': ['Missing data for required field.'], 'memory': ['Missing data for required field.'], 'backend_version': ['Missing data for required field.'], 'max_shots': ['Missing data for required field.'], 'open_pulse': ['Missing data for required field.'], 'n_qubits': ['Missing data for required field.'], 'gates': {0: {'name': ['Missing data for required field.'], 'qasm_def': ['Missing data for required field.'], 'parameters': ['Missing data for required field.']}}, 'backend_name': ['Missing data for required field.']}\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available backends:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<IBMQSimulator('ibmq_qasm_simulator') from IBMQ(ibm-q-ornl, bes-qis, argonne)>,\n",
       " <IBMQBackend('ibmq_20_tokyo') from IBMQ(ibm-q-ornl, bes-qis, argonne)>,\n",
       " <IBMQBackend('ibmq_poughkeepsie') from IBMQ(ibm-q-ornl, bes-qis, argonne)>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit import IBMQ\n",
    "\n",
    "# insert your token & URL here\n",
    "IBMQ.enable_account('<your token>', url='<your url>')\n",
    "\n",
    "# check available backends\n",
    "print(\"Available backends:\")\n",
    "IBMQ.backends()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ibmq_20_tokyo\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "BackendStatus(backend_name='ibmq_20_tokyo', backend_version='1.2.5', operational=True, pending_jobs=0, status_msg='active')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# execute on the IBM Tokyo 20 Qubit Machine\n",
    "backend = IBMQ.get_backend('ibmq_20_tokyo')\n",
    "print(backend)\n",
    "backend.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the UCCSD ansatz circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def get_ucc_ansatz(theta):\n",
    "    circuit = QuantumCircuit(2, 2)\n",
    "    circuit.x(0)\n",
    "    circuit.ry(theta, 1)\n",
    "    circuit.cx(1, 0)\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the naive measurement circuits\n",
    "more automations can be done in here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the number of shots\n",
    "shots = 1000\n",
    "\n",
    "def measure_zi(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    return (counts.get('00', 0) + counts.get('10', 0) - counts.get('11', 0) - counts.get('01', 0))/shots\n",
    "\n",
    "\n",
    "def measure_iz(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    return (counts.get('00', 0) + counts.get('01', 0) - counts.get('10', 0) - counts.get('11', 0))/shots\n",
    "\n",
    "\n",
    "def measure_xx(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.h(0)\n",
    "    circuit.h(1)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    return (counts.get('00', 0) + counts.get('11', 0) - counts.get('01', 0) - counts.get('10', 0))/shots\n",
    "\n",
    "\n",
    "def measure_yy(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.h(0)\n",
    "    circuit.sdg(0)\n",
    "    circuit.h(1)\n",
    "    circuit.sdg(1)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    return (counts.get('00', 0) + counts.get('11', 0) - counts.get('01', 0) - counts.get('10', 0))/shots\n",
    "\n",
    "def measure_hamiltonian(theta):\n",
    "    return 5.9 + .22 * measure_zi(theta) - 6.1 * measure_iz(theta) - 2.14 * measure_xx(theta) - 2.14 * measure_yy(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### run the experiment with different theta value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-3.141592653589793\n",
      "the current values[] array is: \n",
      "[11.77248]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-2.6179938779914944\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-2.0943951023931957\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-1.570796326794897\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-1.0471975511965983\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-0.5235987755982996\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-8.881784197001252e-16\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "0.5235987755982978\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "1.0471975511965965\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996, -0.22271999999999914]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "1.5707963267948948\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996, -0.22271999999999914, 2.59764]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "2.094395102393194\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996, -0.22271999999999914, 2.59764, 6.116479999999999]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "2.617993877991493\n",
      "the current values[] array is: \n",
      "[11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996, -0.22271999999999914, 2.59764, 6.116479999999999, 9.31788]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import qiskit\n",
    "values = []\n",
    "for theta in np.arange(-np.pi, np.pi, np.pi / 6):\n",
    "    values.append(measure_hamiltonian(theta))\n",
    "    # print out the values after each runs in order to save progress from program collapse/network issue\n",
    "    print('theta is: ')\n",
    "    print(theta)\n",
    "    print('the current values[] array is: ')\n",
    "    print(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the results\n",
    "Superimpose the results of the simultaneous measurements and that of the naive measurement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4FNXbxvHvk05I6KFDElqA0EmQIkW6YuEVBRVUUEHsShFREQWsYMOGqICKIoqo/ECp0osQeq8JEAiQAAkJIXWf94/dhNAXSDIp53Ndc5GdnZ25dxL22XNm5oyoKoZhGIbhYnUAwzAMI28wBcEwDMMATEEwDMMwHExBMAzDMABTEAzDMAwHUxAMwzAMwBSEAk9EXhWRb3No3UtE5ImcWLeRc787ESknIstEJF5EPszu9V9nljUi0sfKDMZ5blYHMHKWqr6TG9sRkb7AE6p6a25srzDIwd/dACAGKKY3eSGSiLwKvOp46Aa4A+ccjw+qavDNrN/IXaaFYBiFjz+w40aKgYhc8CVSVd9RVR9V9QEGAqszHptikP+YglBAiMgwETni6AbYLSIdHPPfFJGpjp8DRERFpJ+IHBaR0yIyUERCRWSLiMSKyOdZ1pn52ote73bRtusAE4AWIpIgIrGO+d1EZKOInHFs783LrOtRETkkIjEi8lqW511E5BUR2S8iJ0XkVxEpleX5u0VkuyPzEkeGjOdURGpkeTxFRMY4fi4jIrMdrzslIstF5JL/ByLylYiMu2jeXyIy6Gr7+zLrmSIiX4jIHMey/4lI9SzPf+rYN2dEZL2ItL7c/heRf0Tk2YvWvVlE7nX8XFtEFjje024R6XmlPMCjwMuO31VHEfEUkU9E5Khj+kREPB3LtxORSMf7PQZMvtx6r0ZE2orIBhGJc3QRhV5hucoiskNEnhORh0Vk5UXPvyoi0x0/lxKRn0UkWkTCReRlEZHrzWZcRFXNlM8nIAg4DFR0PA4Aqjt+fhOYmmW+Yv/w9gI6A0nAn0BZoBJwAmh78Wsver2b4/ES7N1EAH2BFRflagfUx/7FowFwHOh+0bq+AYoADYFkoI7j+ReANUBlwBP4GpjmeK4WcBbohL2L4mVgH+DheF6BGllyTAHGOH5+1/H+3R1Ta0Aus0/bOPapOB6XxN4VUvFq+/sy65kCnASaYe9S+Qn4JcvzfYDSjucGA8cAr8v87h4BVmZ5XV0g1rFvijry9HOspzH2LqG6V8k0JsvjUY59XRbwA1YBo7P8DtOA9x3bKnKVv8PL/Q2UBc4APR3Z+gLRQHHH82sc+6AmsB941DG/qON1gVnWtRPo5vj5V+A3wAeoAYQDva3+v5jfJ9NCKBjSsf9nrSsi7qoaoar7r7L8aFVNUtX52D9Yp6nqCVU9AizH/oFy01R1iapuVVWbqm4BpgFtL1rsLVU9p6qbgc3YCwPYux9eU9VIVU3G/uF4n6N10guYo6oLVDUVGIe9qLR0IlYqUAHwV9VUVV2ujk+YiyzHXlgyvrHfh7075CjXv7//UNW1qpqGvSA0yrKPpqrqSVVNU9UPHesNutw6gEYi4u943BuY6dg3dwIRqjrZsZ6NwO/A/U7sj4x1jXL8DUQDbwEPZ3neBoxU1WRVPXfZNVzZPcAmVf3VkW0KEAncnmWZBsAiYKiqfg+gqmcd76EPgIiEYC+c8xytlx7AMFVNUNV9wCcXZTZugCkIBYDjP8SL2D80T4jILyJS8SovOZ7l53OXeeyTHblE5BYRWexo1sdh/5Avc9Fix7L8nJhl2/7AH46unVjs3w7TgXLYv6UfzHiRqtqwf0Ou5ESssdhbE/NF5ICIvHK5hRxF4hfgQcesh7B/mN/I/r7Se0REhojITkd3SixQnEv3EaoaD8wBHnDMejAjD/Z9dUvGvnKspzdQ/iqZsrpgfzp+zvp+olU1ycl1XWvdGevP+rt6FHvr4K+Llvse+/sAe2GY5iiq5bF/dh26yjqNG2AKQgGhqj+r/Qwff+zfbN/PhtWeBbyzPL7aB8zlvmX/DMwCqqhqcexdNc728x4GblfVElkmL0cr5ij29wmAo++4CnDEMSvxSrlVNV5VB6tqNeBuYNCV+v+xt2juc3wrvwX7N9aM9dz0/nYcL3gZe3dKSVUtAcRx5X00DXhQRFpg7/Jb7Jh/GFh60b7yUdWnnIxywf4EqjrmZbiZM5EuXnfG+o9kefwq9q7L7y86nrMU8BKR5tgL4I+O+cewt1qqXmWdxg0wBaEAEJEgEWnvaEonYf+Wb8uGVW8C2ohIVREpDgy/yrLHgcoi4pFlni9wSlWTRKQZ9m/ZzpoAvJ3RRSIifiJyj+O5X4FuItJBRNyx970nY+/7zsj9kIi4ikhXsnRTicidIlLDUUTisLc6LruvHF0vMcC3wDxVzThYnl372xd7/3w04CYibwDFrrL839g/XEcB0x0tI4DZQC3HgVh3xxQqWQ60X8M04HXHPi4DvAFMvcZrnDULaCwi94mIm4g8gv3D+58syyQD/4e9cH+bcXDY0Ur7EZgInFTVMMf8ZOxdaO+ISFHHQfoXsjFzoWUKQsHgCbyH/cPrGPYDeVf78HaKqi4ApgNbgPXYP3iu5F9gO3BMRGIc854GRolIPPYPmV+vY/OfYv8wme94/Rrs39JR1d3YuxA+w/6e7wLuUtUUx2tfcMzL6Dr5M8t6awILgQRgNfClqi7myn4GOjr+zZBd+3seMBfYg73LIwn7t/3LcnwQzrw4j6M7qTP27qSjjkwZB4GdMQYIw/573gpscMy7aap6HHtL7DXsB9efBe5U1biLlktyLFcDmJDljKEfsJ+Y8CMXetLx70Hsf3vfcr4LzbhBGWdQGIZh5Dki4ou99VlbVQ9da3nj5pgWgmEYedlzwBJTDHKHGbrCMIw8yXEhXCL2riQjF5guI8MwDAMwXUaGYRiGQ77qMipTpowGBARYHcMwDCNfWb9+fYyq+l1ruXxVEAICAggLC7M6hmEYRr4iIhdfLX5ZpsvIMAzDAExBMAzDMBxMQTAMwzCAfHYMwTByQ2pqKpGRkSQl3egAn4ZhDS8vLypXroy7u/sNvd4UBMO4SGRkJL6+vgQEBGBuwmXkF6rKyZMniYyMJDAw8IbWYbqMslH58iBy6VTe2VHpjTwhKSmJ0qVLm2Jg5CsiQunSpW+qZWsKQjY6fvz65ht5lykGRn50s3+3piBkMxevFIqEhiEDWoHPsWu/wDAMI48wBSEbqCr/HThJ6Ts3UvmZRZTtMhy3Cqtxv38IiBkrysgeTzzxBDt27MiWdQUEBBATE0NsbCxffvlltqzTyD3vvPNOjqzXFISbcDIhmW+WHaDDR0t5cuICBtSdyjQdRYTrWnZSlHZVf6f0g3/j6mPOVimocvO40bfffkvdunWzdZ2mINyctLQ0S7ZrCkIeYbMpK/fF8OzPG2j37t9snjuJ95PfYUORZ3inyLeUcItkjC2NNJT5eDC+2ovUfWw2i3aaAwkFUU4cNzp79izdunWjYcOG1KtXj+nTpwPQrl27zKFbfHx8GDp0KMHBwXTs2JG1a9fSrl07qlWrxqxZswCYMmUKzz77bOZ677zzTpYsWXLBtl555RX2799Po0aNGDp0KAkJCXTo0IEmTZpQv359/vrLft/7iIgI6tSpQ//+/QkODqZz586cO3cOgP3799O1a1eaNm1K69at2bVrV+Zr2rdvT4MGDejQoQOHDtlvadC3b19mzJiRmcHHxweAqKgo2rRpQ6NGjahXrx7Lly+/ZN8EBAQwfPhwGjVqREhICBs2bKBLly5Ur16dCRMmZC43duxYQkNDadCgASNHjsyc3717d5o2bUpwcDATJ04EID09nb59+1KvXj3q16/Pxx9/fMn+jomJIWMctSlTpnD33XfTvn17OnTocMXtRUREULt2bfr27UutWrXo3bs3CxcupFWrVtSsWZO1a9dm/r4fe+wxmjVrRuPGjTP3+ZQpU7j33nvp2rUrNWvW5OWXX878nZ07d45GjRrRu3fvK/wV3SBVzTdT06ZN1SrHz5zTLxbv1dveX6APD39b/xp5lya9VV51ZDHVcbVV572mtzVdqLzmpbyJeo5ER4301JQ3fPXkiIr64vBhOvLPrZqUmmbZezCcs2PHDqeXhStPN2rGjBn6xBNPZD6OjY1VVdW2bdvqunXrHNtF//77b1VV7d69u3bq1ElTUlJ006ZN2rBhQ1VVnTx5sj7zzDOZ6+nWrZsuXrxYVVX9/f01Ojpaw8PDNTg4OHOZ1NRUjYuLU1XV6OhorV69utpsNg0PD1dXV1fduHGjqqref//9+uOPP6qqavv27XXPnj2qqrpmzRq97bbbVFX1zjvv1ClTpqiq6nfffaf33HOPqqo++uij+ttvv2Vus2jRoqqqOm7cOB0zZoyqqqalpemZM2cu2Tf+/v765Zdfqqrqiy++qPXr19czZ87oiRMntGzZsqqqOm/ePO3fv7/abDZNT0/Xbt266dKlS1VV9eTJk6qqmpiYqMHBwRoTE6NhYWHasWPHzG2cPn36kv0dHR2t/v7+mfu1UqVKmeu60vYy9tmWLVs0PT1dmzRpov369VObzaZ//vln5v4YPnx45r48ffq01qxZUxMSEnTy5MkaGBiosbGxeu7cOa1ataoeOnTogn12OZf7+wXC1InPWHMdwlXYbMryfTFMW3OQE7tXcZesYKbHWkp4nEY9iyF1e0CDnuDfClxcqZ36NCs32khJh2SBN0hmpivM8PLi47NfsSxsBU/tf4FX+9xOjbI+Vr89I4+qX78+gwcPZtiwYdx55520bt36kmU8PDzo2rVr5vKenp64u7tTv359IiIibnjbqsqrr77KsmXLcHFx4ciRIxx3NHcCAwNp1KgRAE2bNiUiIoKEhARWrVrF/fffn7mO5ORkAFavXs3MmTMBePjhhzO/4V5JaGgojz32GKmpqXTv3j1zWxe7++67M993QkICvr6++Pr64unpSWxsLPPnz2f+/Pk0btwYgISEBPbu3UubNm0YP348f/zxBwCHDx9m7969BAUFceDAAZ577jm6detG586dr7mfOnXqRKlSpQCuuL2qVasSGBhI/fr1AQgODqZDhw6IyAW/p/nz5zNr1izGjRsH2E97zmhNdejQgeLFiwNQt25dDh48SJUqVa6Z70bleEEQkUnAncAJVa3nmDcW+03QU4D9QD9Vjc3pLM46FpfEb2GHWfHff7RIXMRwt9X4u0ehLh5IUBeo3xOp2RncvS543erI1aSkp1wwb5Mtmft9fdnQ9hVazh9J6Jmn+eLz+/Dv9jL3NTMXPhmXqlWrFhs2bODvv//m9ddfp0OHDrzxxhsXLOPu7p75t+Pi4oKnp2fmzxn92m5ubthstszXOHN++k8//UR0dDTr16/H3d2dgICAzNdlbAPA1dWVc+fOYbPZKFGiBJs2bXL6/WXNZbPZSEmx/59p06YNy5YtY86cOfTt25dBgwbxyCOPXPL6rO81a6aM966qDB8+nCeffPKC1y1ZsoSFCxeyevVqvL29adeuHUlJSZQsWZLNmzczb948JkyYwK+//sqkSZMuyHnxvitatGjmz1faXkRExCX5Lvd7UlV+//13goKCLnj9f//9d8k+z+ljFrlxDGEK0PWieQuAeqraANgDDM+FHFeVblMW7TzOoO/m883YIbRZ2pPpKc/wgtsfVAmoCXd/jgzdC72mQt27LykGABuf3IiO1EumDQM3QbP+uD23DqneniEuP1N3zj2MnfwLZ5JSLXi3Rl529OhRvL296dOnD0OHDmXDhg03tJ6AgAA2bdqEzWbj8OHDmX3WWfn6+hIfH5/5OC4ujrJly+Lu7s7ixYs5ePDqoyYXK1aMwMBAfvvtN8D+4bZ582YAWrZsyS+//ALYC01GSycgIID169cDMGvWLFJT7f8HDh48SLly5ejfvz9PPPHEDb/vLl26MGnSJBISEgA4cuQIJ06cIC4ujpIlS+Lt7c2uXbtYs2YNYD8+YLPZ6NGjB2PGjMncbtacWY95OLu968n72WefoY67V27cuPGar3F3d8/cb9kpx1sIqrpMRAIumjc/y8M1wH05neNKjsSe44/VOzm9fiZtk5cw1nU7rm42kv3qQeMxSL0eSLGK2bOx4pXweng66dv/IuCvQQw++BS/jZ1D7Yfeo1H1ytmzDSNXlSt3+QPI5crd+Dq3bt3K0KFDcXFxwd3dna+++uqG1tOqVSsCAwOpW7cuderUoUmTJpcsU7p0aVq1akW9evW4/fbbGTZsGHfddRf169cnJCSE2rVrX3M7P/30E0899RRjxowhNTWVBx54gIYNG/LZZ5/Rr18/xo4di5+fH5MnTwagf//+3HPPPTRs2JCuXbtmfttesmQJY8eOxd3dHR8fH3744Ycbet+dO3dm586dtGjRArAftJ46dSpdu3ZlwoQJ1KlTh6CgIJo3bw7YP8D79euX2Rp49913ARgyZAg9e/Zk4sSJdOvW7bq35+rq6lTeESNG8OKLL9KgQQNsNhuBgYHMnj37qq8ZMGAADRo0oEmTJvz0009ObccZuXJPZUdBmJ3RZXTRc/8Dpqvq1Cu8dgAwAKBq1apNr/WNJavy5R3/WX2i4L4HYMZ0SChPufI2fpgfyc5lM6l+7B86uqzHS1JJLFoZz8YP4NqwJ/gFXXP9NyUpjug/X8Vv11QitQwb67/BHfc+gquL6UKy2s6dO6lTp47VMQzjhlzu71dE1qtqyLVea+lBZRF5DUgDrljiVHUiMBEgJCTkuqpX5je3tqOh6gpcOo6ga3p7HqzzO81mrKKznCXRswRpdftA6EN4Vw61n0SeG7yK4/fAF5zd9yCuvz7NXdueZ9X+mdR85DP8KlTNnQyGYRhZWFYQRKQv9oPNHTQnmyk+UdRtPIWHxZ0HG07Hn185a/Mi3r8T6S374F2zA7je2FCx2aFojVvxHrqWbb+NImT315z7ujk7Ql+j7h1P515xMgzDwKIL00SkK/AycLeqJuboxtqO5lkXYQgebFMbDx1sTrn39lP+sam41u5qaTHIIO5e1HvoHaIeXMghtwDqrnuV8A/bk3Jij9XRDMMoRHK8IIjINGA1ECQikSLyOPA54AssEJFNIjLhqiu5UT5R0Ggyb8s5KpDAnS5nmVZpOWc9E3JkczfLv3Zjar68lP9VHUbp+F3wZUtO/vMOpKVc+8WGYRg3KccLgqo+qKoVVNVdVSur6neqWkNVq6hqI8c0MEc23nY0iI0josRkDDIn6dBmdI5sLjt4ebhz12Ovsvme+SylKaX/e5/YT1uihy89ZdAwDCM7FeixjNwCVoPbRd+u3VJwC1xlTaDr0LpJfRoM+otxpd4k8cxJ9LvOpMwaDMnx136xYRjGDSjQQ1ekfn7tCzzysnLFvHjp2Rf5blF7PJe9zcMbviNl9xw87voIat9hdTwjD5kwYQLe3t6XvbLXuH5Tpkyhc+fOVKyYTdcg5RMFuoVQELi6CAM6NaR+/68Z6PkeB+Ld4JcH0ekPQ7z9BjxR8VG0ndKWYwnmhjxWsfp3MHDgwAJZDFT1guE3csuUKVM4evRorm/XaqYg5BNNqpZk3EuP82XQJD5I7UXqzn+IezeUAU0nU+nBUSwLX0GFXqPNPZwtMnrZaFYcWsHopTd/fOpqQ01/8803hIaG0rBhQ3r06EFiov0kvTfffJNx48axa9cumjVrdsG6MgZXW79+PW3btqVp06Z06dKFqKioS7bdt29fnnrqKZo3b061atVYsmQJjz32GHXq1KFv376Zy82fP58WLVrQpEkT7r///sxhG0aNGkVoaCj16tVjwIABmcMxjB8/nrp169KgQQMeeOCBCzJnqFevHhEREURERBAUFMQjjzxCvXr1OHz48BW3dzPDYV9pP8+YMYOwsDB69+5No0aNMvd9oeDMkKh5ZbJy+Ou8wmaz6fS1h7TL65N0xastVEcW06Vv+GiZkaK8VkTxibqpoZeN6xv+WlX16Jmj6jXGPux5kTFFNCo+6qa2f7WhpmNiYjKXe+2113T8+PGqqjpy5EgdO3asqqo2bNhQDxw4oKqq7733no4ePVpTUlK0RYsWeuLECVVV/eWXX7Rfv36XbPvRRx/VXr16ZQ7R7Ovre8HwzRs3btTo6Ght3bq1JiQkZG7jrbfeUtXzw0urqvbp00dnzZqlqqoVKlTQpKQkVT0/vHTWzKqqwcHBGh4eruHh4Soiunr1alXVq27vZobDvtp+zjr0dX5jhr8uRESEnqFVaOJ/P22GVeeeKi/xsezlR4pwh6ShbUbD319YHbNQGb1sNDa1d2ukazqjl47mi2439zu43FDTANu2beP1118nNjaWhIQEunTpcslre/bsyfTp03nllVeYPn0606dPZ/fu3Wzbto1OnTrZc6anU6FChctu+6677socorlcuXIXDN8cERFBZGQkO3bsoFWrVgCkpKRkjuOzePFiPvjgAxITEzl16hTBwcHcddddNGjQgN69e9O9e3e6d+9+zffv7++fOdbQmjVrrrg9uPHhsDOGp77cfi6sTEHIp2qU9SHqj2pMeH47iI2vKMLLrmm833gyLBsBmH6j3BAVH8XkTZMzhz1PSU9h8qbJjGg7gvI+N/47uNxQ02Dv0vnzzz9p2LAhU6ZMueQOaAC9evXi/vvv595770VEqFmzJlu3biU4OJjVq1c7ve0rDS/t6upKp06dmDZt2gWvS0pK4umnnyYsLIwqVarw5ptvZg4bPWfOHJYtW8b//vc/3n77bbZu3XrV4bkvHl76cttzNq86OTx11v1cWJljCPnZre+AKhNIZTqpjMGTVmLL09dZFDRZWwcZMloJOSE+Pp4KFSqQmpp6xVEuq1evjqurK6NHj6ZXr14ABAUFER0dnVkQUlNT2b59+w1laN68OStXrmTfvn2A/RaQe/bsyfxAL1OmDAkJCZlDRmcMv33bbbfx/vvvExcXR0JCAgEBAZlDTW/YsIHw8PDr2p6zbmR46ouHBS8sTAshP6t8/jqL/nqOpvgwzdWdRlUvvRetkTMud1OklPQUVkXmzLUuo0eP5pZbbsHPz49bbrnlih9avXr1YujQoZkfsh4eHsyYMYPnn3+euLg40tLSePHFFwkODr7uDH5+fkyZMoUHH3ww8+5oY8aMoVatWvTv35969epRvnx5QkNDAXv3VJ8+fYiLi0NVef755ylRogQ9evTghx9+IDg4mFtuuYVatWpd9/accSPDU/ft25eBAwdSpEgRVq9eTZEiRZzeP/lZrgx/nV1CQkI046bXRpbhvQFxS6fjo1P5X6VBLItqQ6cJv4OLaQDeCDP8tZGf3czw1+YTIx87duz8Ld1tqa789PEDfOHej04V/yVu8cdWxzMMI58xBaEA8fP15M4n3mCBNqPo8rdJ3H/tA4iGYRgZTEEoYGqVL4b3/RM4qqVJ/PkR0hJOWh0pX8pPXamGkeFm/25NQSiAWtWrzvaWn1As7ST7vnnE3qdkOM3Ly4uTJ0+aomDkK6rKyZMn8fLyuuF1mLOMCqjbu3RjfuSzdD78KWunjabZQ29YHSnfqFy5MpGRkURHR1sdxTCui5eXF5UrV77h15uzjAqw9HQbWz7sRvDZ/9jSeTohrTpZHckwDAuYs4wMXF1dCHryR067lqb8/KfYFX7Y6kiGYeRhpiAUcN7Fy+DWcwrl5RTHfnycE3GF+9J8wzCuzBSEQqB07VZE3zKcdrb/+OPrkSSmpFkdyTCMPCjHC4KITBKREyKyLcu8UiKyQET2Ov4tmdM5CrsKXYcQXeE2+p79jo9/+BWbLf8cOzIMI3fkRgthCtD1onmvAItUtSawyPHYyEki+D08iVSvMvQ5NJJP5piD84ZhXCjHC4KqLgNOXTT7HuB7x8/fA9ceIN24ed6lKNr7eyq7nCRo7atM+++g1YkMw8hDrDqGUE5VM+7fdwwoZ1GOQkeqNocOI+jmupZd//uYFXtjrI5kGEYeYflBZcft3a7YoS0iA0QkTETCzIVC2cO11QukVevIa24/8ulPM9h3ovCN+24YxqWsKgjHRaQCgOPfK96tQlUnqmqIqob4+fnlWsACzcUFtx4TcfHx40M+5pnJSzmZkGx1KsMwLGZVQZgFPOr4+VHgL4tyFF5FS+N2/ySqSDTPnf2CAT+EkZSabnUqwzAslBunnU4DVgNBIhIpIo8D7wGdRGQv0NHx2Mht/i2R217lTpdV1DryO0NnbDEDuhlGIZbjg9up6oNXeKpDTm/bcMKtg+DgSkaHT+XOLTX5uExRBnVy7taEhmHkrKx3RcyqXDn7DbKym+UHlQ2LubjA/03E1bskP/p+wXeLtvDHxkirUxmGweWLwdXm3yxTEAzw8UN6fEuZ1CN8XfInhs3Ywtrwiy8dMQzDCuKZSqn75iL9W4FPDjQLsjAFwbALbI20G86t5/6lv89KBvwYRkTMWatTGUahVyw0nJTa49GKa6DN6BzdlikIxnmtB0NgWwanf0sNPcRjU9YRm5hidSrDKLRcvFLwC13BZy7LCESh8eQcbSWYgmCc5+IKPb7FxbMYPxb7ipOnTzNw6npS0mxWJzOMQqnYLQd4xmccfcWV8riApOdoK8EUBONCPmWhxzcUidvP/6r9yZoDp3jtj63mdFTDyGXR8ck0afY/npFIviWF1ZIObinQeDJlAnKmlWAKgnGpau2g7ctUPfwnX9ffzW/rI/lyyX6rUxlGoTJhyT5GFRnPGYFXOD+SgIdXOj0/z5lWgikIxuW1HQYBrekcMZYBtVMYO283c7ZEXft1hmHctGNxScSvnUorOcMwkjgp51voKekprIpclSPblfzUFRASEqJhYWYc/1wTfwwm3IrNuzS1179LUslkjk9rTkrU+fsZ5dQFMoZRmL09YxUDt/bCp0JNPAcstF8vdBNEZL2qhlxrOdNCMK7MtzzcOxGX6N0Mi51BeoIXZXuE4VosMXORnLpAxjAKq8jTiQRs/pCSkoDnPZ/cdDG4HqYgGFdXvT20HszjDX6i49ZDiKuNsveuBzFnHhlGTvhz9v940GURiY2fgAoNcnXbpiAY19ZuOMsOtuSLdsMpFZbEqSpPU6TpRqtTGUaBE3HiDG1Y5a1EAAAgAElEQVT3vstZj9L4dBmR69s3BcG4Nlc3Hvz9O86levFDy+cRl+24tB8Frma4bMPITmG/f0h9l3Bsnd8Gr2K5vn1TEAynHI2vyMNz36OBdyyf4ElikbkUbb7e6liGUWCER+yn87GJhBcLpXhIL0symIJgOKVcOZjnv4T3bGk8iQf/pwKtx1CuUprV0QyjQIj+fRhepFDivvEgYkkGUxAMp2zcG4VXi8mMcElkLel8Ix6U8pzH0xP/szqaYeR7B9fPo1n8AsIqP0zJqnUty2EKguGU0ctGY1MbaQIPkYgH8APujF/7rrkfs2HcjLQU3OcOIZKyBPccZWkUUxAMp6yOXE1Kun3k0/2iPEMS7XDlKV1vhrUwjJsQNe9DKqYeYkPd4RQvnvsHkrPK8VtoGgXDxicvOs1UFX5/nDe3/UmvNQs5emsgFUsUsSacYeRXsYcoFfYJ/xLKbXc/bHUa00IwbpAIdPsI9a3AR66f8dW8zVYnMox85/TMwaTblCPNR+Lr5W51HGsLgoi8JCLbRWSbiEwTES8r8xjXqUgJ3O77lsoSQ6NtY9h3IsHqRIaRf+yeS8lD8/nO9T56tG9hdRrAwoIgIpWA54EQVa0HuAIPWJXHuEH+LUhqMYgerstZNvNLq9MYRv6QkkjSrMHstVWiaNsX8PbIG733VncZuQFFRMQN8AaOWpzHuAHeHYdzxLcB90V9xM6dW62OYxh5ni7/EK+zkXzkMYCHWtawOk4mywqCqh4BxgGHgCggTlXnX7yciAwQkTARCYuOjs7tmIYzXN0o3nsKIuD6xwBINxerGcYVxexFV37KzPRbadmhO17urlYnymRll1FJ4B4gEKgIFBWRPhcvp6oTVTVEVUP8/PxyO6bhJJ/y1VkXPIJaKTs4/OdbVscxjLxJFf17CInqwWTvx+gZWsXqRBewssuoIxCuqtGqmgrMBFpamMe4SS27D+Rvl3ZU3Po5ejBn7uhkGPna9pnIgSW8n3I/vTuE4umWd1oHcJ0FQURcRCS7rpw4BDQXEW8REaADsDOb1m1YwMvdlXOd3iPSVoak6Y/DuVirIxlG3pF0Bp37Kntdq7Oi+F30aFrZ6kSXuGZBEJGfRaSYiBQFtgE7RGTozW5YVf8DZgAbgK2OLBNvdr2Gte5pFsT7PkNxTzyObfZL9gvYDMOAJe9CwnGGJD7Kcx1r4+5q9Tk9l3ImUV1VPQN0B/7B3uefLZfUqepIVa2tqvVU9WFVNYPi5HNuri5063oXH6X2wGX7TNg8zepIhmG9qC3ofxOY7d6FhDINuadRJasTXZYzBcFdRNyxF4RZjv5+87XPuKLb65VnZbk+bJBgdM4QOGnGOjIKMZsN5gwmxb0Er8Xfy0udauHqYs3w1tfiTEH4GogAigLLRMQfOJOToYz8zcVFGNy1Ls+cG0iKusLvT0B6qtWxDMMam6ZC5Fo+loepWL4Cd9SrYHWiK7pmQVDV8apaSVXvULuDwG25kM3Ix1rXLIN/tZq8bhsARzfA4nesjmQYue/sSVjwBtGlmjIhrhkvdaqFSx5tHYATo52KSAngESDgouWfz6FMRgEgIrzctTb3fnmKvgHdCV7xMVS/DQLbWB3NMHLPojfRpDMMsj1M/Uol6Fy3nNWJrsqZLqO/sReDrcD6LJNhXFWTqiXpVLcc/aLuJb1kNZj5JCSesjqWYeSOw2thww/sDOjD8riyDOpcC7Ho1pjOcqYgeKnqIFWdrKrfZ0w5nswoEIZ0DiI6xY0pFUfA2WiY9Zw5FdUo+NLTYPYg1Lciz0R2pknVErSrlfdHWnCmIPwoIv1FpIKIlMqYcjyZUSAElffl/xpX4oPNXpy59VXYNRvWT7E6lmHkrHXfwPGtLA4cTHi8MLhzUJ5vHYBzBSEFGAus5nx3UVhOhjIKlpc61sKmynun20O122DucIjeY3Usw8gZZ6Lg37dJr96Rl7f707xaKVpWL211Kqc4UxAGAzVUNUBVAx1TtZwOZhQcVUp50/sWf6avP8rBNh+Chzf8/hikmesQjQJo/muQnsKvfs8RczYl37QOwLmCsA9IzOkgRsH2zG018HRzYeyqOLjnCzi2FRaNsjqWYWSv/Yth2+8kt3yJD9am0LpmGUID8k8PuzMF4SywSUS+FpHxGVNOBzMKFj9fTx5rFcjsLVFs82kJof1h9eewb6HV0Qwje6Qlw99DoFQ1JuldnE5MZXDnIKtTXRdnCsKfwNvAKsxpp8ZNGNC2GiW83Rk7bzd0Hg1+deCPpyDB3PjIKABWjoeT+zjb8X2+WnGEjnXK0qhKCatTXRdnrlT+/nJTboQzCpZiXu481bY6S/dEs+ZwIvT4FpLi4K9nzKmoRv52KhyWj4O63fk60p8zSWm81KmW1amumzPDX4eLyIGLp9wIZxQ8j7YMoFwxTz6YuwstF2xvKeydB2u/sTqaYdwYVfjnZXBxI7bNKCatjOCO+uUJrljc6mTXzZkuoxAg1DG1BsYDU3MylFFwebm78kKHWmw4FMuinSeg2QCo2Rnmvw7Ht1sdzzCu3645sHc+tBvOhI3nOJuSxosd81/rAJzrMjqZZTqiqp8A3XIhm1FA3R9SmYDS3oydt5t0Be75EryKw4zHIfWc1fEMw3kpZ+GfYVA2mOi6ffl+VQR3N6xIrXK+Vie7Ic50GTXJMoWIyECcGBTPMK7E3dWFwZ2D2H08nlmbj4CPH3T/CqJ3woI3rI5nGM5b+gGciYRuH/LV8kOkpNt4oUNNq1PdMGe6jD7MMr0LNAV65mQoo+DrVr8CdSsU46MFe0hJs0HNjtD8GVg7EXbPtTqeYVzbiV32U6cb9eFYicZM/e8g9zauRDU/H6uT3TBnuoxuyzJ1UtX+qro7N8IZBZeLizC0axCHT53jl3WH7DM7joRy9eGvpyH+mLUBDeMqos4cZeN37bF5FIVOb/H54r3YbMrz+bh1AM51Gb0gIsXE7lsR2SAinbNj4yJSQkRmiMguEdkpIi2yY71G/tCulh/NAksxftE+ElPSwM0T7vsOUhLhj4H2Ww8aRh4074/HaZx8ll/K1SYyxZvp6w7TM7QKVUp5Wx3tpjjTZfSYqp4BOgOlgYeB97Jp+58Cc1W1NtAQ2JlN6zXyARFhWNcgYhKSmbwywj7TLwi6vgMHFsOaLyzNZxiXc+zkXrqGr+Q/0ul/ZAXvzl2DiPBc+xpWR7tpzhSEjFGZ7gB+UNXtWebdMBEpDrQBvgNQ1RRVjb3Z9Rr5S1P/UnSoXZYJS/cTm5jimNkPat8JC98iet8i2k5py7EE04Vk5A0bfu9LeYTnSSJV0/lxx0c81KwqFYoXsTraTXOmIKwXkfnYC8I8EfEFsqMtHwhEA5NFZKOjO6roxQuJyAARCRORsOhoM8RBQTSkSxAJyWlMWOq43lEE7v4MipYhfUZfNhxcweilo60NaRjA8aObaHd0K9NIZa2kk2pLIcF1Afc1u+SjK19ypiA8DrwChKpqIuAB9MuGbbsBTYCvVLUx9kH0Xrl4IVWdqKohqhri55f37zhkXL86FYpxT8OKTFkVzvEzSfaZ3qU42fVdyibF8Za6M3nTZNNKMCy3f2Y/XIFXScqc5+KifLV+rHWhspEzZxnZVHVDRneO4wK1Ldmw7UggUlX/czyegb1AGIXQS51qkZaujF+0N3PeiIgFfC/pPI0HZW3pppVgWOvYVppHH+AzUoiQ82NvpWsqE/9ZZWGw7ONMCyFHqOox4LCIZIwP2wHYYVUew1r+pYvyYLOqTF93mIiYs0TFRzF502RG6zlcgZfSxbQSDOuowvzXiU0qwdvvR+H+RRz+52ZTYtEueFNJ+2Kj1QmzhWUFweE54CcR2QI0At6xOI9hoefa18DNVfhowR5GLxuNTW2Ei/I9qTyJB+VMK8Gwyr5FcGAJo5YOIzapJCVu3YMtyY0zawvWzSOduQ7hQxEJzomNq+omx/GBBqraXVVP58R2jPyhbDEvHmsVyKzNR1kcvoKUdPtZR2+TnNlKWBVZMJrmRj5iS4cFI6BkIF+uewLPyqfwrnWcM2urYUt2tzpdtnKmhbATmCgi/4nIQMfpooaRI55sU51iXm608PkOHanoSCX8TRvuTR7leddibHxwjtURjcJm009wYgd0fJNUmzslb9tBWrwXZ9YVrNYBOHdQ+VtVbQU8AgQAW0TkZxG5LafDGYVPcW93Brarzr+7TrAu4tT5J1oPAU2HFR9bF84ofJIT4N+3oXIzqHsP3nWO4lkxjthlQWiaq9Xpsp1TxxBExBWo7ZhigM3AIBH5JQezGYVUv5aBlPV13EQn405qJf2hUW9YPwXOHLU0n1GIrP4cEo5Bl7dJSrNRpv1uko8V4+y2ShcsVq6cRfmymTPHED4GdmO/MO0dVW2qqu+r6l1A45wOaBQ+RTxcOTK/JusiTuNdIxoR+7VqAY8MJjXFBss/sjqiURjEH4OVn0Ld7lClGZNWhiM+5/h9RB1UBVUyp2MF5OQ3Z1oIW4CGqvqkqq696LlmOZDJMDi6ogqpp70p2WY3YG8lHIzzZ/Km3rDhe4g7Ym1Ao+Bb/Dakp0LHkcQkJPPl4v10rFOWltXLWJ0sxzhTEDYDQRfdKKe6iLipalxOBzQKKZsLsctr4VHuDN51zncRvb18CKgNVphWgpGDjm+HjVPtt3gtVY1PFu7hXGo6r9xex+pkOcqZgvAlsAaYCHwDrAZ+A3Zn1zDYhnE5iTsrknK8GCVa7wEX+/BZh+KqQuM+sOEHiIu0OKFRYC14Azx9oc0Q9p2IZ9raw/S+pSo1yubfm984w5mCcBRo7LheoCn24wYHgE7ABzkZzijshNNLg3AvmYhv44PnZ7cebO+4NWccGTlh3yLYtxDavAzepXjn7114u7vm61tjOsuZglDLMeQ1AKq6A6itqgdyLpZh2CWF+3EuojTFW+5FPFLtM0uYVoKRQ2zp9tZBCX9o1p+V+2L4d9cJnmlfg9I+nlany3HOFIQdIvKViLR1TF865nkCqTmczyikzp/GJ8QuqYOrdyrFb9l/fn5GK8GccWRkp83T4Pg26Pgm6S4ejJmzk0olitC3ZYDVyXKFMwXhUWAf8KJjOgD0xV4MzMVpRo44duz8KX3Jx4pzd8OKlGsbzqbdjmGHS1SBJg+bVoKRfVLOwr9joFIIBP8fMzdEsjPqDMNur42Xe8G7CO1yrloQHBekfauqH6rq/zmmcaqa6BgWOyGXchqF3NAuQaTblE8W7jk/89ZB9n+Xf2hNKKNgWf0FxEdBl7dJTE1n3PzdNKxSgrsaVLA6Wa65akFQ1XTAX0Q8cimPYVxWlVLe9Gnuz69hh9l7PN4+M7OV8CPEHrY2oJG/xR+HFZ9AnbuhanO+WRbO8TPJjOhWB5GbvmNwvuFMl9EBYKWIjBCRQRlTTgczjIs9174mRT3ceH/urvMzWw+2/2taCcbNWPIOpCdDxzc5cSaJr5ft5/Z65QkJKGV1slzlTEHYD8x2LOubZTKMXFWqqAcD21Vn4c4TrA13DHxXvDI0ecR+EVHsIWsDGvnTiZ32Y1Gh/aF0dT6cv4fUdBuv3F7b6mS5zpnRTt9S1beAsRk/Ox4bRq57rFUg5Yt58e4/O88PfNd6kH2wI3PGkXEjFrwBHr7Q9mV2Rp3h1/WHeaRFAP6li1qdLNc5M7hdCxHZAexyPG7oOPXUMHJdEQ9XXupUk42HYpm7zTGimGklGDdq/2LYOx/aDEGLlOSdv3dSzMud59rXsDqZJZzpMvoE6AKcBFDVzUCbnAxlGFfTo0llapb14YN5u0lNtw9pwa0ZrQRzLMFwki0d5o+wX+jYbABL9kSzfG8Mz3eoSQnvwnkejVP3Q1DVi0/hSM+BLIbhFDdXF4Z1rU14zFl+Wef40yxeCZo8aloJhvO2TIfjW6HDSNJcPHhnzk4CSnvzcHN/q5NZxpmCcFhEWgIqIu4iMgT7bTWzhYi4ishGEZmdXes0Cr4OdcrSLKAUny7cQ0Jymn3mrS+BuMCycdaGM/K+lERYNBoqNYV6PZgedpi9JxJ45fbaeLg59T25QHLmnQ8EngEqAUeARo7H2eUFsrHAGIWDiDD8jtrEJKTwzTLHsFoZrYRNP8Hpg1dfgVG4rfkC4o9C5zEkpKTz8YI9hAaUpEtweauTWcqZs4xiVLW3qpZT1bKq2kdVT2bHxkWkMtAN+DY71mcULo2rluSO+uX5ZvkBTsQ7hrRoPcjeSlhuWgnGFSScsF+EVvtO8G/JhCX7iUlI4bVudQvVRWiX48xZRn4i8qqITBSRSRlTNm3/E+BlwJZN6zMKmaFdapOSZmP8or32GcUqQtO+sOlnOB1hZTQjr1ryLqQlQce3OBp7jm+WH+DuhhVpVKWE1cks50yX0V9AcWAhMCfLdFNE5E7ghKquv8ZyA0QkTETCoqOjb3azRgETWKYoDzaryrS1h9kf7Rha69aXQFzNGUfGpU7sgvXfQ8jjUKYG4+btRoGXuwZZnSxPcKYgeKvqMFX9VVV/z5iyYdutgLtFJAL4BWgvIlMvXkhVJzpuzhPi5+eXDZs1CprnO9TEy82FsXN322eYVoJxJQtHgkdRaDuMrZFxzNx4hMdaBVK5pLfVyfIEZwrCbBG5I7s3rKrDVbWyqgYADwD/qmqf7N6OUfD5+XoyoE115m4/xvqDp+0zM1oJ5owjI8OBpbBnLrQejHqXYsycHZQq6sHTt1W3Olme4UxBeAF7UUgSkTMiEi8iZ3I6mGFcjydaB1LGx5P3Moa0KFYBQvrZWwmnwq2OZ1jNZoP5r0PxKnDLQBbsOM5/4ad4qWNNinm5W50uz3DmLCNfVXVRVS9VLeZ4XCw7Q6jqElW9MzvXaRQuRT3deLFjTdZFnGbhzhP2ma1eBBc3c8aRAVt/hWNboMNIUl08eO+fXVT3sx9/Ms5z5iwjEZE+IjLC8biKiDTL+WiGcX16hVahWpmivD93F2nptiythGmmlVCYpZ6DRaOgYmOo14Of1hzkQMxZXr2jDm6uhfcitMtxZm98CbQAHnI8TgC+yLFEhnGD3F1deLlrEPtOJPDbesdtNW99CVzdzbGEwmzNl3DmCHQeQ1xyOp8u2kvL6qVpX7us1cnyHGcKwi2q+gyQBKCqp4HCOfKTked1CS5Pk6ol+HjBHhJT0sC3PDTtZ795+qkDVsczcltCNCz/GIK6QcCtfLl4H7HnUnmtkN0JzVnOFIRUx72VFewXqmEuJDPyKBHh1TvqcCI+mUkrHN1Et77oaCWY6xIKnaXvQWoidHqLw6cSmbwygh5NKhNcsbjVyfIkZwrCeOAPoKyIvA2sAN7J0VSGcRNCAkrRqW45Jiw9wMmEZHsrIeQxeyvh5H6r4xm5JXoPhE22/+7L1OT9ubtwcYEhnc1FaFfizFlGP2EfXuJdIArorqq/5XQww7gZw7oGkZiSxmf/7rPPaPWCvZVgrl4uPBaOBHdvaPcKGw6dZvaWKAa0rkb54l5WJ8uznL0fwi5V/UJVP1dVMzKpkefVKOtLr9Aq/PTfQQ6ePOtoJTwOm38xrYTCIHw57P4bWr+EepdmzOwd+Pl68mRbcxHa1ZhzrowC68WOtXBzcWHsPMeQFhmtBHPGUcGWcRFasUrQ/Gn+3nqMDYdiGdypFkU93axOl6eZgmAUWOWKefFE60Bmb4li8+FY8C1nbyVsmW5aCQXZthkQtQk6vEGyePDe3J3ULu/L/SFVrE6W55mCYBRoA9pUo1RR+5WpqupoJXjAsrFWRzNyQsZFaOUbQP2e/LDqIIdPnePVO+rg6mJOM70WUxCMAs3Xy53n29dg9YGTLNkTbW8lhJpWQoH13wSIOwydx3D6XBqf/buXtrX8aFPLjJTsDFMQjALvoVv88S/tzfv/7CLdltFK8DSthAIkKj6Ku79rhW3ZOKjVFaq15dNFe0lITuO1bnWsjpdvmIJgFHgebi4M6RzErmPxzNwQCT5lz7cSYvZZHc/IBqOXjabT4Q1oylnoNIrwmLNMXXOQXqFVqVXO1+p4+YYpCEah0K1+BRpWLs5HC/aQlJpuWgkFSFR8FMs3TmYg7kxySeNYkeK8989OPN1ceKlTTavj5SumIBiFgouL8MrtdYiKS2LKqgh7K6HZE/ZhkWP2Wh3PuAljlo7i03RXkoBRksqz/3uNeduPM7Btdcr6movQrocpCEah0aJ6aW4L8rMPcJaYAi1NKyG/i4qPosiGH2ivrgwiiUhbCn/s+YlSvok80bqa1fHyHVMQjEJl2O21iU9O44vF+8DHz9FK+M20EvKpb+YNYbTNlf+RyrekAmDTdMpWmk0RD1eL0+U/piAYhUrt8sXo0aQy3686yOFTifZWgpsXLP3A6mjG9UpP5b7d8zmL0p8kyLjMQNI4lrTF0mj5lSkIRqEzqFMtROCjBXvsrYTQJ+xXt0bvsTqacT2WfkDd1BQGTP+J42/ZKDZ3L/7nZuP5XQzH3txodbp8yRQEo9CpWKII/VoF8uemI2w/Ggctn7e3EpaZVkK+ERkGyz/k+00P8seuu3HxTqZ48/0k7i1H8uHSHD9udcD8ybKC4Lg382IR2SEi20XkBauyGIXPU+2qU7yIO+/9s8txLKE/bJ3BiYiVtJ3SlmMJx6yOaFxJylmYOQCKVeT5ue8DSqnO2xC3dE4vqW11unzNyhZCGjBYVesCzYFnRKSuhXmMQqR4EXeeva0Gy/fGsGJvjL2V4O5NxKynWXFoBaOXjrY6onEl80fYb4fa/SvOJBenRNvdFA06xumltUk75WN1unzNsoKgqlGqusHxczywE6hkVR6j8Hm4hT+VShTh3X92YitSmoRGDxFyKpz6Npi8abJpJeRFexdC2HfQ4hkIbI1Pg0MUb76f+I1ViV8XaHW6fC9PHEMQkQCgMfDfZZ4bICJhIhIWHR2d29GMAszTzZWhXYLYfvQMszYf5a200xwDplEEL1u6aSXkNYmn4K9nwK8OtB/Byn0xlOq8jXMH/Di1IJjzpxkZN8rygiAiPsDvwIuqeubi51V1oqqGqGqIn58ZsdDIXnc3rEhwxWK8M28Vn2/7md4kEoQLH6W7mFZCXqIKcwZB4km492v2nkpl4NT16Bkfov9qDHrhR1m5chblzOcsLQgi4o69GPykqjOtzGIUTvYhLWqz8+xkUtPTWSLpjCKFvnjwYLqaVkJesXUGbP8D2r1CjG9t+k1Zh6ebK6veC8GW7I4qF0zHTB2/IVaeZSTAd8BOVf3IqhyG0bqmH8lp+0hX+5Wuo0lmMWmMt7mzdP4Sa8MZEHcE/h4MlZuRdMtz9P8hjJiEZL57NITKJb2tTlegWNlCaAU8DLQXkU2O6Q4L8xiF2aR1+J+bTYl/d2J7S+n94T4Sz5bh55Ti9rtwGdaw2eCvpyE9DVv3CQyesZ1Nh2P5pFdjGlYpYXW6AsfKs4xWqKqoagNVbeSY/rYqj1G4pZ4oTsL2iviGhONWPJGohAo8/MdEGpTbDnOHWx2v8Fr3DRxYAl3eZlxYKnO2RjH89tp0rVfe6mQFkuUHlQ0jr4hdFoSmu+DXYx3ikcq8/R15f+ULsH4ybPvd6niFT/RuWPAG1OzCr9qBL5fs58FmVelvRjHNMaYgGIZD+hlvYv5sinups/jdsxHExuv/joDKoTDrBfvFUEbuSE+1X43s7s26Bm/x6h/baF2zDKPuCcZ++NHICaYgGEYWSQfLcGpBPYpUi6Zkhx2k2dzhvkng4gK/9YO0ZKsjFg7LxkLUJqLavMvjvx+iml9RvujdBHdX85GVk8zeNQwuPG89YXNV4tYGUqzpQSq2jYASVeGeLyFqEyx806qIhUfkelg2jqS699NzeVk83Fz57tFQinm5W52swDMFwTCwn7ee9Tz2mIV16FinHJ4ttrNk9wmocyfcMhDWfAm7zLkPOSYlEf4YgPqW54kTPYmOT+bbR0OoUsqcXpobTEEwjMtwdRE+faARtcsX49mfN7L7WDx0GgUVGsKfT0HsYasjFkwL3oCT+/is2CBWRKbycc9GNDKnl+YaUxAM4wqKerrxXd8QvD1ceWzKOmKSgPsmgy0dfn/cfuDTyD77FsG6bwir8AAf7avAK7fX5vb6FaxOVaiYgmAYV1GheBG+fTSEk2eTGfBDGEnFAuCuT+Dwf7D4bavjFRyOgevifKrTO7wrD4RW4ck25vTS3GYKgmFcQ4PKJfikVyM2HIrl5Rlb0Ho9oGlfWPEx7FtodbyCYc5gbAnRPHL6cUJrVGR093rm9FILmIJgGE7oWq8CL3cNYtbmo3y6aC90fQ/K1oWZT8KZKKvj5W9bZ8D2mXyh95FYup45vdRCZq8bhpOealudHk0q88nCvfy1/ZT9eEJqIszsbz+uYFy/uCPYZg9im0sQU13/j0l9QylexJxeahVTEAzDSSLCu/fWp1lgKYbO2ML6c+XgjnEQsRyWjbM6Xv5js5H+59OkJCfxUspTTHj0FnN6qcVMQTCM6+Dh5sKEPk2pUNyLAT+Ecbhqd2jwACx9D8KXWx0vX9F13+AavoTRqb15qVdXGlctaXWkQs8UBMO4TqWKejCpbyip6TYe/yGM+I7vQalq8PsTcDbG6nj5Q8xe0uaNYHF6Qyp1fJo7zOmleYIpCIZxA6r7+fBVn6YciD7LszP2knbvJDh3Gv540j6Gv3Fl6amcmtqXhHR3VtZ9k6fa1bA6keFgCoJh3KBWNcowuns9lu6JZsx6N+j6rv001FXjrY6Wpx3+axSlYrfxQ+kXGdaznTm9NA8xBcEwboJ9fP5ApqyK4IfU9lC3O//f3n2HV1XneRx/f5IQIAFCMZhAAJG2MKFJUBkVsKzAjDso6zqow4iMorI4uiiWwY5lfLCX0UelDMrawFVhUWRtoIDSlSKIKAISDCChRFJuvvvHPUhQSiTlJJfv63l4kntyzrmfH0nu956S749374L1n4YdrUrasGw26Z89wcwavRl8+fmpr3gAAA4XSURBVLV+e2kV498N58ropn7RRnh3TF3B7Pa3QUoGTB4SPYXkfrJt+3aKpwxlCw1oP+Rpv720CvKC4FwZ7W2E1y6tHsMmr2HdmU/Czmx4Y3i0daojvyjCJ88Mp7l9x44+j5GR7heRq6JQC4KkvpJWSVoj6aYwszhXFsk1Exh7SRa1E+O5eHohu3reCl9Mg0+fCTta6MyMsf8cR7+8qaxtNYi2Pc4JO5I7iNAKgqR44EmgH9ABuFBSh7DyOFdWTepHG+Ft2ZXPn5d3I9K6D7xzC3y3OOxooXrq7YWc9+29bEtqyfEDx4Qdxx1CmEcIJwJrzGytmRUALwH9Q8zjXJl1yqjPwxd0YdH6XEYxDEtOjU69uWdH2NEqTVoaSNF/dTI30HTOrRxjuQycOB5q1A47njuEMAtCU6DkLCMbgmX7kTRU0gJJC3JyciotnHNHql/HdEb2acdLy3YzucWdsP1bmHrNUXM9YfPm6MfENss5b0BP+sfP4a5ZNzJzeddwg7nDqvIXlc3sGTPLMrOs1NTUsOM4VyrDekcb4Y2cn8SK9n+F5a/Bon+GHauSGEntN9L83/+Lx+K2MG93I+6bdV3YoVwpJIT43BuBZiUeZwTLnKv2JHHvgEzWb8vjvM+6M7/5adR760bI6A7H/ibseBXCzHh35fekX7qKhMZf80TiHBKJY1DNbCLJW2BXWtgR3WGEeYQwH2gjqaWkRGAg8GaIeZwrVzUT4nl6UDfSUpI4f/MlRBLrRa8nFOwOO1q5m/PVFgY8NYfLJi6gW+0VvMZt9FE817OHNYpAz9FhR3SlEFpBMLMiYDgwA1gJvGJmy8PK41xFaJicyNhLupMdqcfN+iu2ZTVMvyHsWOVm8bc/cPFz87jo2U/I357Ne61fZWrqjaTX/IELyeNpCiGhALqOhzrZYcd1hxHmKSPMbDowPcwMzlW01o2jjfAuGRfhlNSL6L/kBWh5GnQeGHa0I/ZF9g4efGc1M1dsJjUpnhc7LeXkdU+h7/J4aEt7bm+wmF0JRfs2UITafUcTvdPcVVVV/qKyc7FgbyO8EZv78k2dLjBtBN+v+5heE3qRvav6vHNet3U31760mH6PzmbeV1t58MRdzGt0Jz1W34+angBXzeX5jvnsSijYf8OEAtqdNSec0K7UZNXoVrisrCxbsGBB2DGcO2J3T1vB1I8W8mHdW9iSAO32bGBw1pU8+fuq/c45O3cPj733Ja/MX09CvBieVYeh+RNIXDEZUppBn3uh/b9F//jAVTmSFppZ1mHX84LgXOWJFBvH/Wkhp7eYzMSaTzGJQq4vguxH13JschrZVexgYdvuAp76YA0T566j2Iw/ZaUzIuV96s57ECIFcMo1cOoISPSpL6uy0haEUK8hOHe0iY8TGyd3YcrVd9OqZiG3U4M/xhvTh5zB2BljIHI2xIffBXTnnkKem/01Yz/6mryCIs7rmsGNbTfR+KOhsHQ1tOkTnf+hUauwo7py5EcIzlUy1d0E1xwPNfbQ2uIYQg0GWyLpEpGkVOK6XIROGATHtKn0bHsKI0yc+w3/+OArtucV0i8zjRtOTqblonthxRvQ4Djoez+061vp2dyR81NGzlVROmcYdB0bvR0zEG8JDIh04zKSOCN+MQkUs6FuZ7b/y4U0PvkCGjdqVKGZCoqKeWXBeh5/70s278inZ9tURp7Zgo7rJsKsB6MrnXYd/PZqqFGrQrO48uenjJyrqjLm7lcMACIq4tUf8ug/6GXuW7uW9HWvc8aOGWTOv4mdn97J6/GnsCLtXOq1OomOzRrQqWkKDZITyxwlUmy8uXQjD8/8km+35ZHVogGPDuzKyUUL4Y0rYNtaaP8H6HMP1G9e5udzVZsfIThXyQ51I07JX8cf84tYt/Rd4ha/QIvsd6hpe1hVnMErkd68FjmVug3T6JiRQueMFDo2rU/HjBTq1Pzle7y0tH0N50o8E+ndN5N54SpWb95Fh/R6jOzTjt6Nd6O3b4bVb0GjNtDvfmh9ZrmM24XHTxk5V0Ud+AUajj2Wg99ltGcHLJtCZOFE4jctIqIElib/lhfye/H6znYUE4cErVLr0KlpCp0yUujUrD4d0utROzG+xI6MWi22Ur/nF9Rsksvxqclc96/t6NeuHnFzHoWPHoG4BOh9I5x0FSSU/SjEhc8LgnOxavMKWPw8LH0JftxGpG5T1jc/l/drn83HW5NZuiGXnJ35ACTEibxNdcnflEJhTj2S2mZTq8VWinJrs/3jNmxf1ISEL6fD23+D3G8h83w4ezTUaxLyIF158oLgXKwryodV02HR8/DVe4BBy15Y10FkNz2LpZvy+Xzjdh4Yl0tiWi7xtQuJ7E4kd05rdi5tTtv6a1n18A3RbRt3gN+NgeNODXtUrgJ4QXDuaLJ9PSz5b1jyQnRCnlr1odMF0HUQatIJMOLS11Lc91KSXx/PLSdMZESPJ0lMqg2nj4Lul0G832MSq7wgOHc0Ki6Grz+MnlJaORUiBSz8rjNjFw/ixSbz6NPlBR4obEhGYh7jF1/MpZPugDqNw07tKpgXBOeOdnnb4PNXWfb882Qe8zkRM+IlFpox/MWX+XpHvyrXKsNVjNIWBO926lysSmoIJ11B5uOzua9DXx6Li3A5P3JqQj4nPDzNi4H7BS8IzsW4TbuyuevLNxlBHs+pkD3FBYxfMr5atd12lcMLgnMxbvSs0RRb8X7LIhZh9Ic+raXbnxcE52Lc3A1zKYjs3yqjIFLAnA0+YY3bn99n5lyMW3zF4rAjuGrCjxCcc84BIRUESWMkfSHpM0n/I6l+GDmcc87tE9YRwkwg08w6AauBm0PK4ZxzLhBKQTCzd8ysKHg4D8gII4dzzrl9qsI1hCHAWwf7oqShkhZIWpCTk1OJsZxz7uhSYa0rJP0fkHaAL40yszeCdUYBWcAAK0UQSTnAunINus8xwJYK2ncYYm084GOqLnxMVU8LM0s93Eqh9TKSNBi4AjjTzPJCCVGCpAWl6fVRXcTaeMDHVF34mKqvUP4OQVJf4AagV1UoBs4558K7hvAEUBeYKWmJpKdDyuGccy4QyhGCmbUO43kP45mwA5SzWBsP+JiqCx9TNVWt5kNwzjlXcarCbafOOeeqAC8IzjnnAC8IP5E0OuittETSO5KahJ2prGKxZ5Sk/5C0XFKxpGp9G6CkvpJWSVoj6aaw85SVpHGSvpe0LOws5UFSM0nvS1oR/MxdE3amiuYFYZ8xZtbJzLoA04Dbwg5UDmKxZ9QyYAAwK+wgZSEpHngS6Ad0AC6U1CHcVGU2AegbdohyVARcZ2YdgJOB/4yB79EheUEImNmOEg+TgWp/tT0We0aZ2UozWxV2jnJwIrDGzNaaWQHwEtA/5ExlYmazgG1h5ygvZrbJzBYFn+8EVgJNw01VsXyCnBIk3QP8GcgFTg85TnkbArwcdgj3k6bA+hKPNwAnhZTFHYak44CuwCfhJqlYR1VBOFx/JTMbBYySdDMwHLi9UgMegV/RM6oImFSZ2Y5UacbkXGWRVAeYAlz7szMJMeeoKghmdlYpV50ETKcaFITDjSnoGXUO0Z5R1eI02K/4PlVnG4FmJR5nBMtcFSKpBtFiMMnMXgs7T0XzawgBSW1KPOwPfBFWlvJSomfUH7xnVJUzH2gjqaWkRGAg8GbImVwJkgSMBVaa2UNh56kM/pfKAUlTgHZAMdEW21eaWbV+xyZpDVAT2BosmmdmV4YYqcwknQc8DqQC24ElZtYn3FRHRtLvgEeAeGCcmd0TcqQykfQi0Jtoq+jNwO1mNjbUUGUg6VRgNvA50dcFgL+Z2fTwUlUsLwjOOecAP2XknHMu4AXBOecc4AXBOedcwAuCc845wAuCc865gBcEFzMk1Zc0rMTj3pKm/cp9DK7oTreS7pB0fUU+h3NHwguCiyX1gWGHXevQBgPVvvW5c0fCC4KLJX8HWgVzWowJltWRNDmYF2JS8NenSOom6UNJCyXNkJQu6XwgC5gU7KO2pNskzZe0TNIze7ffS1KKpHWS4oLHyZLWS6oh6fJg26WSpkhK+nlgSR/snddB0jGSvgk+jw/ms5gfzGdxRbA8XdKsIN8ySadV0P+lOwp5QXCx5CbgKzPrYmYjg2VdgWuJzjlwPHBK0J/mceB8M+sGjAPuMbPJwALg4mAfPwJPmFl3M8sEahPtC/UTM8sFlgC9gkXnADPMrBB4Ldi2M9HWyX/5FWP5C5BrZt2B7sDlkloCFwX77wJ0Dp7buXJxVDW3c0elT81sA4CkJcBxRFteZAIzgzf88cCmg2x/uqQbgCSgIbAcmPqzdV4G/gi8T7Qn0T+C5ZmS7iZ6KqsOMONX5D4b6BQctQCkAG2I9kAaFxS1183MC4IrN14QXKzLL/F5hOjPvIDlZtbjUBtKqkX0xT3LzNZLugOodYBV3wTuldQQ6Aa8FyyfAJxrZkuDrrO9D7BtEfuO1EvuW8DVZvaLIiKpJ/B7YIKkh8xs4qHG4Vxp+SkjF0t2AnVLsd4qIFVSD4i2OJb0mwPsY+8L9JagJ/75HICZ7SL6zv1RYJqZRYIv1QU2Be/mLz5Ilm+IFhF+tv8ZwFXBtkhqG1yfaAFsNrNngeeAE0oxXudKxY8QXMwws62SPg4meX8L+N+DrFcQnIp5TFIK0d+DR4ieDpoAPC3pR6AH8CzReZyzib7oH8zLwKvsfxRwK9EZtnKCjwcqVg8Ar0ga+rO8zxE9vbUouJCdA5wb7H+kpEJgF9EZ/pwrF97t1DnnHOCnjJxzzgW8IDjnnAO8IDjnnAt4QXDOOQd4QXDOORfwguCccw7wguCccy7w/5YreMtDX2trAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# simulated results copied from the other ipython notebook\n",
    "values_simulated = [12.318439999999999, 13.70312, 12.840320000000002, 10.260440000000001, 6.42872, 2.7654000000000014, -0.4627999999999991, -1.6761599999999992, -1.1323599999999996, 1.9101599999999999, 5.18572, 9.35652]\n",
    "# physical results from Tokyo 20 for naive measurement method\n",
    "values_physical = [11.77248, 12.465000000000002, 11.226320000000001, 8.612680000000001, 5.894, 2.0026800000000016, -0.6282799999999995, -1.3829999999999996, -0.22271999999999914, 2.59764, 6.116479999999999, 9.31788]\n",
    "# array of theta values\n",
    "thetas = [-3.141592653589793, -2.6179938779914944, -2.0943951023931957, -1.570796326794897, -1.0471975511965983, -0.5235987755982996, -8.881784197001252e-16, 0.5235987755982978, 1.0471975511965965, 1.5707963267948948, 2.094395102393194, 2.617993877991493]\n",
    "# physical results from Tokyo 20 for simultaneous measurement method\n",
    "values2_physical = [11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001, 2.16244, 5.65808, 9.46308]\n",
    "\n",
    "# plot them!\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(thetas, values2_physical, 'bs', label = 'simultaneous measurement')\n",
    "plt.plot(thetas, values2_physical)\n",
    "plt.plot(thetas, values_physical, 'g^', label = \"naive measurement\")\n",
    "plt.plot(thetas, values_physical)\n",
    "plt.legend(loc='best')\n",
    "plt.title('simultaneous vs naive for Tokyo')\n",
    "plt.xlabel('theta values')\n",
    "plt.ylabel('energy sums')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the difference of simultaneous vs naive for Tokyo 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd8W/W5+PHPI+8VjzjbtpyQBEjSLGeyW2gZpUDLCJRCuqDctr+2t+W2dNLblI7L7bi0tAVaQqDQhE3KKC0tJUCchEzIHo4cO3GG5dhxPOKh7++PcxSEsGPZlnQ0nvfrpZeto6OjR9KRHn23GGNQSimlQuFyOgCllFLxQ5OGUkqpkGnSUEopFTJNGkoppUKmSUMppVTINGkopZQKWdImDRF5SER+bP9/rojsCLjtdBHZKCLNIvIVEckSkb+KSJOIPOFc1AMnIheISK3TcThBRI6LyLgIHDdpX9P+EpFyETEikhqBY78kIgsjcNyzRWSXff5cFe7j9yOOTPu1K3EqhkBJmzQCGWNeN8acHrDpm8Crxpg8Y8w9wDXACGCoMeZaR4J0ULx/ORpjco0xVZF+HBHxiMhFkX4c9V7GmEuNMUsicOgfAb+1z59nB3MgO7Edty+dItIRcP0PYYo3KsKe9ROEG1gadH2nMaarvwcSkdSB3E8p5Tg3sGUgdwz+3BtjLg247SGg1hjzvUFH6ARjTFJcgBnAeqAZWIaVFH5s33YB1psI8C+gG2gHjgN/ATqATvv65+z9PgtsA44CLwPugMcywJeAXcBee9sZwD+ABmAHcF3A/g8B9wIv2PGtBk4LuH1ywH0PAd+xt7uAO4A9gBd4HCjq5flfANQC3wHqAQ9wY8DtGcD/Avvsx/gDkAXkAG2Az37+x4HR9rZi+77fBbqAIfb1RcCvT3XcgMe9HNgINAIrgakBt3mA24G3gSb7fcvs5fmNB16z96sHlgW9H+MDXuvfAS/Zz+VNYCTwa/u93A7M6Om+Affv6bx5xH6N2uzjftPe/gRw0I5rBTC5H+/7qc6ZfOBh4AhQDXwPcNm3/RD4c8C+5fbzSLWvfxqosh9zb+B5EPSazgEq7femDvgtkB702tyGdZ432s9F7NtS7Pe93n6sLwXG0MNj9fpeA4XA8/ZzPWr/XxJw338Dn8c61xqBKQG3DbPfk+F9nW9B8ewJej8zsM775fb7sRu4JWD/HwJPAn8GjgGfP8V30UPY51DQ9i/x7mf5aWCEvT3Tfu1K7OsfBGqAs4E/AXcFHefvwH/Y/38AeN1+vm8Dlw76uzRcX8qxfAHS7Q/WfwJpWNVNnfTw4Q88CYNOiMAP4ZX2SXMmVmnte8DKoA/TP4Ai3v3irQE+Y+8/w/4wTQo4ibxYH9JU4FFgqX1bHtYH9hv2yZMHzLVv+yqwCiixT+r7gL/08hpcgPXF/kt73/OBFuB0+/Zf2R+IIvsx/gr8tKfXx962Arg64CTd4z8h7ds+HsJxZwCHgblYXzILsb48MuzbPcAarA9rEVaSvq2X5/cXrOTlsl+nc4Lej8CkUQ9U2Pv9C+uL82Y7hh9jVU32K2kExHtRUFyftZ93BlZi2hh0rN7e977OmYeB5+xjlwM7efcHzQ/pJWnYxz0W8L6PIiCRBcVeAcyz71duv/5fC3ptngcKgDKsL/VL7Ntuw0rApfZ79yp9J40e32tgKHA1kG0/3yeAZ3v6vAIPEvAlivVF/LdQzrdeYroo4PoKrB8cmcB0+/l+KOA17wSuwjoHs3o6ZvA5FLDtMqwfF1Pt498P/N2+7WTSAK7A+gE2w77tPKzz15+sRwOt9muYae/7DazvvYuxEuDYQX2fhvPLOVYv9gt7wP/C2ttWMvCk8RL2B9S+7rLfKHfAh+lDAbcvAF4Piuk+4M6Ak+iPQSfQdvv/G4ANvTyvbcCFAddH2Sfu+z6YvJs0cgK2PQ58HxCsBBL4K3c+75aS3vP62NsWAfdgfaEcxEpgP7NP1DasD3pfx/09sCjouDuA8+3/PcCnAm77H+APvbwWD9sftJIebgtOGg8E3Pb/gG0B1z8ANPZ034D7h5w0guIosI+XH8L73us5g/WF14GdQOzbvgD8u5fztZz3Jo1GrC/hXr/Yeon/a8AzQa9NYHJ+HLjD/v9fBCR44CP0nTRCfa+nA0d7+rwCFwF7Am57E7g5lPOtl5gusv8vxaqByAu4/afAQwGv+YoQX8eT51DAtkeBHwWdKz6sUrA/aXwLK0GcEbCfYJXkzrWv3w48bf//Yawfy4Hfe8/436OBXpKlIXw0sN/Yr5qtehDHcwP/JyKNItKIVVwVYEzAPjVB+8/172/f50asE8LvYMD/rUCu/X8p1q/43uJ4JuCY27BO7BG97H/UGNMScL0a67UZhvUrbl3Asf5mb+/Na1hfmjOBd7BKVudj/TLdbYzxhnBcN/CNoNel1I7Jr7fXJdg3sd6DNSKyRUQ+e4rYDwX839bD9d4eo19EJEVEfiYie0TkGNaXEEBxwG69Pb9TnTPFWL8cA8/hat57/vXIfv8XYJUE6kTkBRE5o5f4J4rI8yJy0I7/J0Gxnyr+0bz3MxDK563HY4lItojcJyLVdhwrgAIRSenhGK8C2SIyV0TKsRLMM/ZtoZxvvRkNNBhjmoOeU2+f+f4aTcBrZIxpxCoRBh7/61g/BrYH7GewfjB9yt70KayqUv8x9/XwvdfneXIqyZI06oAxIiIB28oGcbwa4AvGmIKAS5YxZmXAPiZo/9eC9s81xvxHiI/VW3fRGqwqocDjZhpj9veyf6GI5ARcL8MqgdVjfVlODjhOvjHG/wVggg+EVVI7Hfi4/dy22se7DCuhEMJxa7CqEgLjzzbG/KWvFyWYMeagMeYWY8xorF/dvxOR8f09Tg9asRKf38jeduT9r9MnsaoyL8Jqgyi3twt9O9U5U49VonQH7F8G+N/3llPFbIx52RjzYayS6XbggV5i+L19+wRjzBCs9rBQYgfrM1caFN9AfQPrXJtrx3Gevf19sRhjurFKPDfYl+cDvugHc74dAIpEJC9gW+BrDj1/TkJ1gID3U0QKgCFBx/848CkRuS3ovg8D14hIBdZr/kLAMYNf9+CY+y1ZkkYlVtXMV0QkTUQ+gVWPPFB/AL4tIpMBRCRfRE7VFfd5YKKI3GQ/fpqIzBaRM0N4rOeBUSLyNRHJEJE8EZkbEMddIuK24xgmIlf2cbz/FpF0ETkXq1HwCWOMD+uL41ciMtw+1hgRudi+zyFgqIjk+w9ijGkF1mHVGfuTxEqsX7Cv2fv0ddwHgNvsX4UiIjki8tGgD2ZIROTagH7sR7E+wL7+HqcHG4FP2qWGS7BKU705xHsTfB5wAqvdIhvrl3qoej1nAr4Y77LPBzf2r9CAmM8TkTL7Pfu2/6AiMkJErrR/PJzAquPu7XXKw/q1e9wujYTyI8fvcazPW4mIFGJ12BioPKwfH40iUoRVRXcqj2GVpm60//cb8PlmjKnBOr9/Kta4ianA53j3NR+svwC3iMgUEcnEqur9lzEmsPS1D7gQ+E5gSdpY3cm3AouxOoB02De9Drjs745UEfkwVjXh44MJNCmShv0ifgKr10gD1gn19CCO9wzwc2CpXVzeDFx6iv2bsd6s67Gy/0H7/hkhPFYzVt3kx+z77cLqPQHwf1iNzH8XkWasRvG5PR3HdhDrC/UAVh3qbQFF3W9hNe6vsp/TK1i/7rD3+QtQZRfr/cX517CqSdYEXM/Dqj4ghOOuBW7B6pVz1N7v0329Jr2YDawWkeNYr8lXTXjGZnwV67X3Vw+dqr/+T4Hv2a/R7Vi/AKuxftltxXp/QhLCOfP/sEoUVcAbWF+OD9r3/QdW76O3sRL78wGHdmElmANYn4Xz6T0Z3I5VWmrG+sJdFmr89v4vA5uwei0O+POG1YEgC6uEtQqrirNXxpjVWK/NaKz2R//2wZ5vN2CVFg9gVXndaYx5pR/3P1XMz2OdP8vt448EbuphvyqsxPEjEQm8fQlWe9wjAfu2Y/0wvAbrh8svgQWD/Vz4W9yVUkrFKRH5CPA7Y0w4qmRPKSlKGkoplahEJB34ClbvwYjTpKGUUnFKRKZjVbXlYQ2ujPxjavWUUkqpUGlJQymlVMgSbsLC4uJiU15e7nQYSikVV9atW1dvjDnVgF4gAZNGeXk5a9eudToMpZSKKyIS0iwZWj2llFIqZJo0lFJKhUyThlJKqZBp0lBKKRUyTRpKKaVCpklDxb265jrOf+h8Dh4/2PfOSqlB0aSh4t6iFYt4Y98bLHptkdOhKJXwNGmouFbXXMdDG54gq/MCFm9YrKUNpSJMk4aKa4tWLCKz4wKKO79OatdULW0oFWGaNFTcqmuuY/HGxUi3tSR6VsfHWLxRSxtKRZImDRW3Fq1YhM/4SDWjAMjyTUe63FraUCqCNGmouFVZW0lHdwdpZhStrjX4aCW743JW1q50OjSlEpYmDRW3NnxhA+3f7SKNEXz7gzdy6zmTGWI+yPMLNGkoFSmaNFRcq2lowxgoL87mM+eMBWDxG3sdjkqpxKVJQ8W1am8LAO6hOYwpyOLyqaP4y5p9NLV1OhyZUolJk4aKax5vKwBjh+YAcMu542jp6Gbpmn1OhqVUwtKkoeKap76FIZmpFGSnATBlTD5nnTaUxW966OjyORydUolHk4aKax5vC+XFOYjIyW23njeOg8fa+eumAw5GplRi0qSh4lq1txW3XTXld/7EYZw+Io8HXq/CGONQZEolJk0aKm51dPmoPdpK+dDs92wXET5/7li2H2zm9V31DkWnVGLSpKHi1v7GNnyG95U0AK6YPprheRk88HqVA5Eplbg0aai45bG72waXNAAyUlP49NnlvL6rnq0HjkU7NKUSliYNFbeq698do9GTG+e4yU5P4Y9a2lAqbDRpqLjl8baSk55CcW56j7fnZ6exYHYpyzcdoK6pLcrRKZWYNGmouFXdQ3fbYJ89eywGeOhNT9TiUiqRadJQccvjbaW8l6opv9KibC77wCgeW72P5nadWkSpwdKkoeJSV7ePmoZW3D00gge75dyxNJ/oYumamihEplRi06Sh4tKBxna6fKbPkgbA1JIC5o0r4sE399LZrVOLKDUYmjRUXPKcnN2275IGWFOL1DW188LbdZEMS6mEp0lDxSX/lOjlxX2XNAAumDic8cNzuX+FTi2i1GBo0lBxyeNtJTPNxfC8jJD2d7mEW84dy9a6Y6zc441wdEolLk0aKi5Ve1soH3rq7rbBrpw+huLcDO5foYP9lBooR5OGiFwiIjtEZLeI3NHLPteJyFYR2SIij0U7RhWbQuluGywzLYXPnF3OazuPsP2gTi2i1EA4ljREJAW4F7gUmATcICKTgvaZAHwbONsYMxn4WtQDVTGn22fY523FXRxaI3igG+eWkZWWwgMrdB1xpQbCyZLGHGC3MabKGNMBLAWuDNrnFuBeY8xRAGPM4SjHqGJQXVMbHd2+fpc0AAqy0+2pRfZzsKk9AtEpldicTBpjgMDRVrX2tkATgYki8qaIrBKRS3o6kIjcKiJrRWTtkSNHIhSuihXV9rrgoXa3DfbZs8fS7TM8tNITxqiUSg6x3hCeCkwALgBuAB4QkYLgnYwx9xtjZhljZg0bNizKIapoe3dK9P6XNADKhmZz6ZRRPLq6muMnusIZmlIJz8mksR8oDbheYm8LVAssN8Z0GmP2AjuxkohKYtXeVtJTXYwckjngY3z+3LE0t3ex7C2dWkSp/nAyabwFTBCRsSKSDlwPLA/a51msUgYiUoxVXaX9JZOcp74Fd1E2Llfo3W2DzSgrZE55EQ++sZcunVpEqZA5ljSMMV3Al4GXgW3A48aYLSLyIxG5wt7tZcArIluBV4H/MsboyKwkV+1t7XXhpf649bxx7G9s48XNB8MQlVLJIdXJBzfGvAi8GLTtBwH/G+Dr9kUpfD5DdUML500sHvSxPnTGcMYNy+H+FXv42NRR/RooqFSyivWGcKXe41BzO+2dvrCUNKypRcaxef8xKqu0AKtUKDRpqLjiqbe62w6051Swj88YQ3FuOg/o1CJKhUSThoor1f2cEr0vmWkp3Dy/nFd3HGHnoeawHFOpRKZJQ8UVj7eVtBRhdEFW2I75qXluMtNc/PF1LW0o1RdNGiquVHtbKC3KJmUQ3W2DFeWkc21FKc9uOMDhYzq1iFKnoklDxZWBzG4bis+dM5ZOn48llZ6wH1upRKJJQ8UNYwzV3pawtWcEKi/O4ZLJI/nzqn206NQiSvVKk4aKG0eOn6C1o5uxIS7x2l+3nDeOprZOHl+rU4so1RtNGipuvDu7bWSSxsyyQma5C/mTTi2iVK80aai4sbfeP7tt+Kun/G45bxy1R9v42xadWkSpnmjSUHGj2ttCqksYE8butsEuOnMEY4tzeGBFFdYsNkqpQJo0VNzweFspKcwiNSVyp22KS/jcOWPZVNvEmr0NEXscpeKVJg0VN6yeU5Fpzwh09cwSinLSeUAH+yn1Ppo0VFwwxlBd3xrR9gy/rPQUbprn5pVth9l9+HjEH0+peKJJQ8WFhpYOmk90RaWkAXDzfDcZqTq1iFLBNGmouOCxu9tGaoxGsKG5GVxTUcLT6/dzuFmnFlHKT5OGigvhnt02FP6pRR6prI7aYyoV6zRpqLjgqW/BJVBSGL2kMW5YLh8+cwSPrKqmtUOnFlEKNGmoOOHxtjKmMIv01OiesreeN47G1k6eXFcb1cdVKlZp0lBxodrbEpHZbftS4S5kRlkBf3x9L90+HeynlCYNFRc83taotmf4iQi3njuOfQ2t/F2nFlH9UNdcx/kPnc/B44l13mjSUDGvsbWDprZOR0oaAB+ZPBL30Gzu06lFVD8sWrGIN/a9waLXFjkdSlhp0lAxzxPh2W37kuISPn/OWDbWNLK2+qgjMajYZ4zh8LF23thVz6/+uYknV2VQeOJrLN7wSEKVNlKdDkCpvvi7244tjn71lN81FaX88h87uX9FFbPLixyLQ8UG7/ET7DjUzK5Dx9lp/91xqJmmts6T+2RyNi6yONH9Lxa9toh7P3qvgxGHjyYNFfM89a1IlLvbBvNPLfKbV3ez58hxThuW61gsKnoaWzvYeTIxNJ9MFN6WjpP7DMlMZeKIPD46dRQTh+cydEgn1z97Hu3dbZS2P4arczyLNy7m++d/n5G5Ix18NuGhSUPFPI+3hdH5WWSmpTgax03zy/nDiip+868tbGj7JsuuWZYQXwLJpK65juufuv59792x9k52HTr+nsSw81Azh5tPnNwnNyOV8cNzuejMEUwYkcvEEXmcPjKP4XkZiMjJ/b74whfpkqP4pINOqSHDdybtZnnClDY0aaiY54nQuuD9NSwvg6tnjmHZWg81Ge8kzJdAMvnhq3exZu9Bbl22mLkjrzxZiqhreneqmKy0FCaMyOXcCcOYOCKXiSPzmDgij9H5me9JDr2prK2ko9sqiZxwbSerew4dXR2srF0ZsecVTZo0VMyr9rZy8eTY+EV/xYw8/rImhZzOSxOqyiEZ7K6v5aWV5zCCS3l7F2yv8jB+eB5zxxYxYUQep4+wkkNJYRYuV9/JoTcbvrDh5P9L1+zjjqffYc+XmhmXIFWamjRUTGtq66ShpSMqU6KH4pGtd9OeMpy8rss47HtWSxtx5I6XHsDFHI6m/onO9PUsnHk5v7/8txF9zFnlhQCsqz6aMElDu9yqmLbP4e62geqa61i8cTFNKc+RQgFpHXNYvHFxQnWnTFR1zXW8trMOQyfNqS/SZqpZsunBiL9344pzyc9KY/2+xOmqrUlDxTSP3d223MHutn6LVizCZ3y0uzbRKTXkdX2MbtOdcIO3EtGiFYtI657MCdcOjFiN29F471wuYWZZAWs9mjSUioqTU6IXOV/SONnAKYbm1BfIMKdDZ3nCNHAmsjerN5DmG0u76+2T2zq6o9M4XeEuZNfh4zS1dva9cxxwNGmIyCUiskNEdovIHafY72oRMSIyK5rxKed5vK2MHJJJVrqz3W3BauA0dxrMnYb9332anPQUvjR52XsaPlVs+sUFzyKk8NJn7jn5Hpo7TVTeu5luq11jfU1ilDYcSxoikgLcC1wKTAJuEJFJPeyXB3wVWB3dCFUs8NTHRnfbYHmZaVxdUcLzm+qoP36i7zsoR1VWeUlPdTGjrCDqjz2tpIAUl7A+QaagcbKkMQfYbYypMsZ0AEuBK3vYbxHwc0DX3ExCHm+rYxMV9uXm+W46un0se6vG6VBUHyr3eKkoK3RkgGhORipnjspjXTImDRFxiciQMD32GCDw01Zrbwt8vJlAqTHmhT7iulVE1orI2iNHjoQpPOW04ye6qD9+AncMNIL3ZPzwPM4ZX8yfV1XT1e1zOhzVi8bWDrYdPMb804Y6FsMsdxEbaxoT4jzpM2mIyGMiMkREcoDNwFYR+a9IByYiLuCXwDf62tcYc78xZpYxZtawYcMiHZqKEn8jeKyWNMAqbdQ1tfPKtkNOh6J6saqqAWNwNGnMdBfS2tHN9oPNjsUQLqGUNCYZY44BVwEvAWOBm8Lw2PuB0oDrJfY2vzxgCvBvEfEA84Dl2hiePKpPjtGIzZIGwIVnjmBMQRYPrfQ4HYrqxaoqL5lpLqaVRL89w6/C/e4gv3gXStJIE5E0rKSx3BjTCYRjJZq3gAkiMlZE0oHrgeX+G40xTcaYYmNMuTGmHFgFXGGMWRuGx1ZxwD9GIxYG9vUmxSV8ap6bVVUN7EiAX5GJqHKPl1nuoqivLx9odH4mI4dkJk3SuA/wADnAChFxA8cG+8DGmC7gy8DLwDbgcWPMFhH5kYhcMdjjq/hXXd9KcW4GuRmxPdvN9bNLyUh18XClx+lQVBD/uhdOVk2BtWxwhbswOZKGMeYeY8wYY8xlxlINfDAcD26MedEYM9EYc5ox5i572w+MMct72PcCLWUkF4+3xdGFl0JVmJPOFdNG8/T6/e9ZhEc5b/XeBgDmjXM2aYDVrrG/sY2DTfHdETSUhvACEfmKiPxSRO4RkXuwGqiViqhqb2tMV00FWnhWOW2d3Ty5rtbpUFSAyj1estNTmFqS73QozPIP8ovzeahCqZ56ESgH3gHWBVyUipi2jm4OHmuPmdlt+zJlTD4V7kIeqfTg84WjyU+FQ2WVl9nlRaSlOD9j0qTRQ8hMc8V9FVUolcWZxpivRzwSpQJUN8R+I3iwm+e7+erSjazYdYQLTh/udDhJ73BzO7sPH+eaihKnQwEgLcXF1JIC1sZ50ggl/T4iIreIyCgRKfJfIh6ZSmqeequ7bSyP0Qh26ZRRFOdmsES738aEVVVWe8b8GGjP8KtwF7JlfxPtnd1OhzJgoSSNDuBuoJJ3q6a0QVpFlH9gX1mcVE8BpKe6+OTcMv698wie+hanw0l6lXu85GWkMnl0uCaxGLyKskK6fIa3a5ucDmXAQkka3wDG2+MlxtqXcZEOTCU3j7eVopx08rPSnA6lX26cW0aKCH9eVe10KElvVZWX2WOLSI2B9gy/mQkwyC+UV3M30BrpQJQKVO2Nzdlt+zJiSCaXTBnJ42traO3ocjqcpHWwqZ299S0xVTUFUJSTzrjinIRPGi3ARhG5z9/l1u52q1TEVHtbGRtH7RmBFp5VzrH2Lp7dcMDpUJLWqiov4Ox8U72pcBeyft9RjInPXnahJI1ngbuAlWiXWxUF7Z3dHGhqi6ueU4FmuQs5c9QQHq70xO0XQ7yr3ONlSGYqZ46KnfYMvwp3IQ0tHeyN03avPrvcGmOWRCMQpfxqGloxJjbWBR8IEeHTZ7n51lPvsGZvA3NjrIokGVRWeZk7bigpLnE6lPcJnLxw3LBch6Ppv1BGhO8VkargSzSCU8nJc3J22/gsaQBcMW0M+VlpLKn0OB1K0tnf2Ma+htaYa8/wO21YLkMyU+N2ZHgog/sCpyLPBK4FdJyGiph319GIz5IGQFZ6Cgtml/KnN/ZS19TGqPwsp0NKGpV7Yrc9A8DlEmbG8eSFoUxY6A247DfG/Br4aBRiU0nK420hPyuNgux0p0MZlJvmufEZw2Or9zkdSlKp3OOlMDuN00fkOR1KryrKCtl56HhcTnAZSvXUzIDLLBG5jdBKKEoNSLW3Na5LGX6lRdlceMZw/rJmHye64ncEcDwxxrCqysvcsUNxxWB7hp+/XWNDHFZRhdJ76hcBl58CFcB1kQxKJTePtyWu2zMC3Ty/nPrjHbz4Tp3ToSSF2qNt7G9si9mqKb9ppQWkuIT1cVhFFUrvqbCsnaFUKDq6fOw/2sbHZ8TGJHODdc74YsYV57BkZXXCPKdYFuvtGX45GamcOSovLicvDKV66qsiMkQsfxSR9SLykWgEp5JP7dFWfCa+G8EDuVzCzfPdbKxpZFNNo9PhJLzKKi/FuelMGB77XVkrygrZWNNIV7fP6VD6JZTqqc8aY44BHwGGAjcBP4toVCppxcO64P11dUUJOekpPFyp81FFkjGGyj3W+AyR2G3P8JvpLqS1o5vtcba2fChJw//qXwY8bIzZErBNqbB6d0r0xChpAORlpvGJmSX89e0DeI+fcDqchOXxtnLwWHvMjs8IVhGnK/mFkjTWicjfsZLGyyKSB8RXeUrFjWpvC3kZqRTlxHd322ALz3LT0eVj6Vs1ToeSsOKlPcNvTEEWI4ZkxN14jVCSxueAO4DZxphWIB34TESjUknL423FXZwdF9UL/TF+eB5njx/Ko6uq464OO15UVnkZlpfBuOL4qNoUESricJBfKIP7fMaY9caYRvu61xjzduRDU8moOoG62wa7eX45B5raeWXbIadDSTj+8Rnz46Q9w6/CXUTt0TYOHWt3OpSQxc7qJCrpdXb7qD3allDtGYEuPGM4YwqyWLJSG8TDbc+RFo40n4ibqim/ijhclEmThooZBxrb6PKZuFoXvD9SU1x8ap6byiovOw/FV4+ZWFfpXz8jThrB/SaNGkJGqiuxkoaI/EJEJkcjGJXc/LPblsdJnfRALJhdSnqqi4crPU6HklBW7fEyKj8z7lZ7TE91Ma2kILGSBrANuF9EVovIbSKSH+mgVHLy1PvHaMTXB78/inLSuWLaaJ5ev59j7fE3WV0siteuCVoSAAAgAElEQVT2DL+Z7kK2HGiivTM+5icLpSH8j8aYs4GbgXLgbRF5TER0ehEVVh5vC9npKQzLzXA6lIhaOL+c1o5unlxb63QoCWHnoeN4WzqYF2ftGX4V7kI6uw3v7G9yOpSQhNSmISIpwBn2pR7YBHxdRJZGMDaVZKq9rbiH5sTlr8X++EBJPjPLCnhkVTU+ny4HO1iVe+qB+GvP8Iu3xvBQ2jR+BezAGtz3E2NMhTHm58aYjwEzIh2gSh4eb0vC9pwKtvCscvbWt7Bi1xGnQ4l7q6oaGFOQRWlRfJ47RTnpjCvOYa0nQZIG8DYwzRjzBWPMmqDb5kQgJpWEun2GmobWhB2jEezSKaMozs3Q+agGyeczrNrrjbuutsFmugtZv+8oxsR+yTOUpLEJOD1oMabTRCTVGBMflXAq5h1obKOz2yRNSSM91cUn55bx6o7DJ5e3Vf23/WAzja2dcVs15VfhLqShpeNkD8JYFkrS+B2wCrgfeACoBJ4Adgx2inQRuUREdojIbhG5o4fbvy4iW0XkbRH5p4i4B/N4KnZVJ0F322A3zi0jRYQ/r9LSxkCdHJ8R5yWNeGrXCCVpHABmGGNmGWMqsNoxqoAPA/8z0Ae2G9fvBS4FJgE3iMikoN02ALOMMVOBJwfzeCq2+adET9SBfT0ZMSSTi6eMZNlbNbR2dDkdTlyq3OPFPTSb0QVZTocyKOOH5ZKXmZowSWOiPR06AMaYrcAZxpiqQT72HGC3MabKGNMBLAWuDNzBGPOqPUkiWKUdXfosQVV7W8hMczE8L7G72wZbOL+cY+1dPLfxgNOhxJ1un2H1Xm/cV02BtVjXzLLCuFj+NZSksVVEfi8i59uX39nbMoDBjE4aAwTOE11rb+vN54CXerpBRG4VkbUisvbIEe2NEo/21rfiLsrB5Urs7rbBZpcXcuaoISxZ6YmLRtBYsvXAMZrbu+K+aspvlruQnYebaWqL7UGfoSSNhcBu4Gv2pQr4NFbCiMoAPxH5FDALuLun240x99vVZ7OGDRsWjZBUmFmz2yZHI3ggEWHhfDfbDzbzVpx0uYwVlVXW+Ix5CVDSAKtdwxjYEOOLMp0yadjtDn80xvzCGPNx+/K/xphWe8r044N47P1AacD1EntbcAwXAd8FrjDG6LJnCcjnM1Q3tCZVI3igK6ePYUhmKktWepwOJa6sqmpgXHEOI4ZkOh1KWEwrLcAlxHwV1SmThjGmG3CLSCSWUXsLmCAiY+3jXw8sD9xBRGYA92EljMMRiEHFgIPH2uno8iVlSQMgKz2FBbNL+duWgxxsip91FZzU1e1jzd6GuJ06pCc5GamcOWoI6+K5pGGrAt4Uke/bXWC/LiJfH+wDG2O6gC8DL2NNivi4MWaLiPxIRK6wd7sbyAWeEJGNIrK8l8OpOJaMPaeC3TSvHJ8xPLZau9+GYvOBYxw/0ZUQjeCBKtyFbNzXGNOrO6aGsM8e++IC8sL54MaYF4EXg7b9IOD/i8L5eCo2JeMYjWBlQ7P50OnDeWzNPr70ofFkpKY4HVJM868HnijtGX4V7kIerqxmx6FmJo+OzQnF+0waxpj/BhCR7IDur0qFjcfbQnqqi1EJUjc9UDefVc7CB9fw0jsHuWrGqToSqsoqLxOG5zIswbpozyyzBvmtrz4as0kjlAkL54vIVmC7fX2a3e1WqbCorm+lrCg76brbBjt3fDHjinNYUulxOpSY1tntY62nIWG62gYqKcxixJAM1sZwY3gobRq/Bi4GvADGmE3AeZEMSiWXZJrd9lRcLuGm+W427Gvk7dpGp8OJWW/XNtLa0Z1w7RlgdcGucBfG9MjwkNbTMMbUBG2KjyWmVMwzxuDxtiTN7LZ9ubqihOz0FJas1Abx3qyqagBgbgImDbCqqGqPtnHoWGz2pAsladSIyFmAEZE0Ebkdq7eTUoN2uPkE7Z0+LWnYhmSm8YmZY/jr2wfwHtdhST2p3OPljJF5FOVEYiSA8/yTF8bqeI1QksZtwJewpvjYD0y3rys1aO+uC64lDb+F88vp6PKxbG1wAV+d6OpmbXVDwvWaCjR5dD7pqa6YraIKZY3wemPMjcaYEcaY4caYTxljvNEITiW+k91tNWmcNGFEHmedNpQ/V1bHdH99J2yqaaK905eQjeB+6akuppXkx+wgv1B6Tw0Tke+IyP0i8qD/Eo3gVOLzeFtIdQmjC5K7u22wm+eXc6CpnVe26UQIgSr3eBGBeWMTN2mAtZLf5v1NtHfGXvNxKNVTzwH5wCvACwEXpQat2mt1t01NCalPRtK46MzhjM7P5OFKj9OhxJTKqnomjRpCfnaa06FE1Cx3EZ3dhnf2x97iqKF8UrONMd8yxjxujHnKf4l4ZCopeJJ0dtu+pKa4+NR8Nyv3eNl1qNnpcGJCe2c36/c1JmRX22AzywqA2FzJL5Sk8byIXBbxSFTSMcbgqdfutr1ZMKuU9FQXv39tK+c/dD4Hjx90OiRHrd93lI6uxG7P8Buam8HY4py4TRpfxUoc7SJyTESaReRYpANTia/+eActHd3a3bYXQ3Mz+NjU0Ty38SBvVq9n0WuLnA7JUauqGnAJzB5b5HQoUeFfyS/WFucKpfdUnjHGZYzJNMYMsa8PiUZwKrFV27PbupN4osK+XD49l25fKtldF7B44+KkLm2s2uNlyph8hmQmdnuGX4W7EG9Lx8kehrEilN5TIiKfEpHv29dLRWRO5ENTic6j3W379MSuX9Dh2k1u18V0+7qTtrTR1tHNhpqjSdGe4ecf5Bdr81CFUj31O2A+8En7+nHg3ohFpCKirrku5urFq70tpLiEMQVZTocSk+qa61i8cTHNKS+TbsZCV1nSljbWVR+ls9sk1KJLfZkwPJe8zNSYa9cIJWnMNcZ8CWgHMMYcBRJz/H4CW7RiEW/seyOmfql6vK2MKcgiPVW72/Zk0YpF+IyPlpQV+Ggnt+sjdJvkLG1UVtWT4hJmlydHewZYE1j62zViSSif1k57rXAD1mA/QIepxpEDxw7wxKpu8jo+GVO/VKu9LUm98FJfKmsr6ejuwEgLrSkryek+n84uYWXtSqdDi7rKPV6mluSTmxHKunGJo8JdyM7DzTS1dTodykmhJI17gGeA4SJyF/AG8JOIRqXC6pZlj5DTeQX5XddiuvNj4peqMYa99Tol+qls+MIGzJ0Gc6fh+c99Bxc5PPGx3Wz4wganQ4uqlhNdvF3blFTtGX4V7kKMgY01sTNVfii9px4Fvgn8FKgDrjLGPBHpwFR4/GvHbjbvOZ121xZAyOi4KCZKG42tnTS3d+kYjRDNHVtE+dBslr2VfJMYrq0+SpfPJMX4jGDTSgtwSWwN8gt1PY3txph7jTG/NcbotOhxoqm1ky89toFuaeRI+o9pc60lr+tiun3ieGljr93dVksaoRERrptdyhpPA1VHjjsdTlRV7vGSliInexMlk9yMVM4YOSSm2jW0BTJBGWO4/clNtJ1I50j6T/FJM82pz5NCIWkdcxyvFz85RkNLGiG7ZmYJKS7h8bW1TocSVZVVXqaVFJCdnlztGX4V7kI27DsaMzMea9JIUH98fS//2HqIH1w+lRP/vQNzp6Hlh2sZV5zDh0f92PF6cU99KyJQWqTdbUM1fEgmHzx9OE+uq6UzRr5AIq25vZPN+5uSsmrKb1Z5IS0d3eyIkTnINGkkoLWeBn72t+1cOmUknzm7/OT2WFqDutrbwuj8LDJSUxyNI94smF1K/fETvLo9OaZMf8vTQLfPJGUjuN/MsthayU+TRoLxHj/Blx/bQElhFj+/Zioi8p7bY2UNao+3lfJibc/orw+ePozheRlJ0yBeucdLeoqLmUnYnuFXUpjF8LyMmGkM16SRQHw+w38+vomG1g7u/eTMHufoiZU1qKu9LTp9yACkpri4uqKEV3cc5tCxdqfDibjKKi8zygrITEveEqmI1QkgVlby06SRQO59dTcrdh7hhx+bzJQx+b3ud7PDa1A3tXZytLVTk8YAXTerFJ+BJ9cldoN4U1snWw4cS+r2DL8KdyE1DW0cjoEfCpo0EsTK3fX86pWdfHzGGG6YU3rKfSfaa1A/umqfIz0yqhv8Pae0emogxhbnMHdsEY+vrcHni61ps8Npzd4GjCGp2zP8/NVzsVBFpUkjARw+1s5Xlm5k3LBcfnzVlPe1Y/Tk5vnl7G9s458ONKjurbfHaOgUIgN2/ZxSqr2trN7b4HQoEVO5x0tGqovp9ip2yWzy6CGkp7o0aajB6+r28f/+soGWE138/saZ5IQ4N49/DeolKz2RDbAH/vUByoq0pDFQl04ZRV5mKsve2ud0KBFTWeWlwl2oPeyAjNQUppXkx0S7hiaNOPerV3ayem8Dd318ChNG5IV8v9QUFzfOc2YNao+3hVH5mUnduDlYmWkpXDl9NC9tPhhTk9mFy9GWDrbVHdOqqQAz3YVs3t9Ee2e3o3Fo0ohjr+44zL2v7uGGOaV8YmZJv+9//WxrDeqHK6Pb/bba26rtGWFw/ewyTnT5WL5xv9OhhN3qvV4AbQQPUFFWSGe3YfP+JkfjcDRpiMglIrJDRHaLyB093J4hIsvs21eLSHn0o4xN+xvb+M9lGzlz1BDu/NjkAR3Dvwb1U+trOdYevV+r2t02PKaMyWfSqCEsTcAxG5V7vGSlpTC1RNsz/GKlMdyxpGGv0XEvcCkwCbhBRCYF7fY54KgxZjzwK+Dn0Y0yNnV0+fjyY+vp6jb87saZg6rmWXiWm9aObp6OUvfN5vZO6o93aCN4mFw/p5QtB445/usz3CqrvMwqL9QFugIU52ZQPjQ7eZMGMAfYbYypMsZ0AEuBK4P2uRJYYv//JHChhNI1KMH9/G/b2bCvkf+5ZipjB/nlO7WkgOmlBTxcWR2V7pvVJ9cF1+qpcLhy2hjSU10JNUK8/vgJdh46rlVTPZjpLmRd9VGMca6rtZNJYwwQeKbX2tt63McY0wU0ARE7k2JxHe1gf9t8kD+9sZdPn1XOZR8YFZZjLjzLTVV9C2/srg/L8U7FnzR0dtvwyM9O49IpI3l2437HG0jDZXWV1Y1YG8Hfr8JdiLel4+TnyAkJUfYTkVtFZK2IrD1y5MiAj/Oj12JvHe1A1d4W/uuJTUwrLeA7l50ZtuNe9oFRFOem83ClJ2zH7I3HqwP7wm3B7FKa27v42+bY/bHTH5VV9eSkp5xyVoNkNcttrZHuZBWVk0ljPxA4dLnE3tbjPiKSCuQD3uADGWPuN8bMMsbMGjZs2ICCqfLWsnzlJLI7L2bxhiUxV9po7+zmi4+ux+US7v3kjLDW9WakpnDDnDL+uf0wNQ2R/QXjqW9heF5G0q6NEAnzxg6lrCibpQkyZqNyj5fZY4tIS0mI37RhNWF4LnkZqY6O13DyXXkLmCAiY0UkHbgeWB60z3Jgof3/NcC/TIQq8xa99kt8tDC080sUtv6U//zrbyLxMAO26PmtbDlwjF9eN42SwvD/Sv/k3DJcIvx5VWS731Z7W7XnVJi5XMKC2aWsqmrAY4+2j1eHj7Wz50iLVk31wuUSZrgLHZ0m3bGkYbdRfBl4GdgGPG6M2SIiPxKRK+zd/gQMFZHdwNeB93XLDYe65jqWbvs9B9Pv4EjazxHfECo3ncV/PFoZExOEPbdxP4+u3sdt55/GhWeOiMhjjMrP4uLJI1j6Vg1tHZGrG/d4W7RqKgKunlmCS+BxhyahDJfKKh2f0ZeKskJ2HGqOajf5QI6W/4wxLxpjJhpjTjPG3GVv+4ExZrn9f7sx5lpjzHhjzBxjTFUk4li0YhE+4wOB1tTXOZB5G8fTnuLlzUf40C9e44EVVY6tlLb7cDPffvod5pQXcftHJkb0sW6eX05TWyfLN0VmsFhrRxeHm09od9sIGJn/7qp+sbIs6ECsqvKSl5nK5NHantGbCnchxsCGfc4spKaVhkBlbSUd3R0nrxtpx5u6mNwx9zC7vJC7XtzGJb9eweu7Bt7IPhCtHV188dH1ZKWlcM8NM0iNcB3v3LFFnD4ijyUrqyPSpe/d7raaNCLhutmlHG4+wWs7o3uehtOqqgbmji0ixZX0Pet7Nb2sAJc41xiurZHQ53rZ/9x2iB89v5Wb/rSGSyaP5HuXnxmRdoVAxhi+9+xmdh0+ziOfncvI/MyIPh5Yi70sPKuc7zzzDuuqjzKrvCisx6/WnlMR9aEzhlOcm8HSt2oiVo0ZSQeb2tlb38KNc8ucDiWm5WakcsbIIY61a2hJIwQXnjmCl792Hrd/ZCKv7TzChb94jf97ZVdE+8U/vraGp9fv56sXTuCcCcURe5xgV80YTV5mKksiMB+V5+QYDU0akZCW4uLqijH8a/thDjc73xbXX5VV1jihedoI3qcKdyEb9h2l24H1VDRphCgzLYUvf2gC//zG+Vw0aQS/emUnF/3yNf6+5WDYq3K2HjjGD57bwjnji/l/H5oQ1mP3JTs9letmlfLSO3Vh7wTgqW+hODedvB6WoVXhcd2sUrp9hqfWxd8khpV7vORnpTFp1BCnQ4l5Fe5CWjq62XEwujNUgyaNfhtdkMW9n5zJY7fMJTs9hVsfWcfCxW+x58jxsBy/ub2TLz22noLsNH59/XRH6nZvmuemy2d4dHV4+/1bPae0PSOSThuWy5zyIp5YW+PoVBMDUVnlZe7YIlzantGnCv/khQ6M19CkMUBnnVbMC185lx9cPokN1Ue55Ncr+OlL2zh+omvAxzTGcMdT77CvoZXf3DCT4tyMMEYcuvLiHC44fRiPrdlHR1f4euLolOjRcd3sUqrqW3jL4/yCPaGqPdpKTUObdrUNUUlhFsPyMhxp19CkMQhpKS4+e85Y/nX7BVw1fQz3vVbFh/733zy7Yf+AfuU9XFnNC+/U8V8Xn86cseFthO6vhWeVc6T5BH/bEp6R8e2d3dQ1tWvPqSi47AMjyc1IjasR4pV7dHxGf4gIFWWFrK2O/nK/mjTCYFheBndfO41nvngWo/Iz+dqyjVx3XyVbDoQ+XfWmmkZ+/MJWLjxjOLeeOy6C0Ybm/AnDcA/N5uEwLQe7r0EbwaMlOz2VK6aP5sV36hwbANZflVVeinLSmTg89NUnk12Fu5CahraoD0DWpBFGM8oKeeaLZ/Pzqz/AniMtfOw3b/D9ZzfT2Npxyvs1tXbyxUfXMzwvk19cNy0m6nRdLuGmeW7WVh8Ny1oN/uktBjuVuwrNglmltHf6WL7xgNOh9MkYw+qqBuaN0/aM/qgot9o11ke5XUOTRphZ8wCV8eo3LuDm+eU8urqaD/7vv3l0dXWP3eOMMXzjiY0cbm7n3htnUpCd7kDUPbu2opSstJSwzH57ckr0Ik0a0TC1JJ8zRubFxbQiNQ1t7G9s0/mm+mny6CGkp7qiPshPk0aE5Gen8cMrJvPCV85lwog8vvvMZq68942Tb7B/7Y5fvLKJV7Yd5ruXncn00tha2jI/O42rZozhuY0HONpy6tJSXzzeFgqz08jP1u620SBiTWL4dm0TWw8cczqcU9LxGQOTkZrC1DH5mjQSzZmjhrDs1nncc8MM6ps7uPr3K/nG45v43iv/w5q9Xu79Zw0f/cAoFp5V7nSoPVp4lpsTXb5B/2K1ek5pKSOaPj7DWtUv1ksblXu8FOdmMH54rtOhxJ0KdyGb9x+L6gJcmjSiQES4Ytpo/vmN8/mPC07juY37+cfqsyg+8W065RBfv2QEsbqK7RkjhzB3bBGPrOq5ei1Ue+tbdInXKCvITufiySN5ZkPsrupnjKGyysu8cUUx+xmIZTPdhXR0+6K6RrwmjSjKyUjlW5ecwdkz/0VHyjaEdI5m3s0vV/3E6dBOaeFZ5dQebePV7YcHdP8TXd0caGrTkoYDFswqpamtk5fD1HU63NZUezh07ASTx2i15UDMLLMH+UWxikqTRpTVNdfx+I57OZR+JzWZN9DKThZvXBxzKwUG+vCkEYwcksmSSs+A7l/T0IYxUF6sJY1oO+u0oZQUZsVsFdWiV5YBsOrwIw5HEp+G5WXgHpqtSSORnVy7A0CsKoNu0x2z65KDNYjxxrllvL6rnt2H+z9dyruz22pJI9pcLuG6WaW8udsb8aV8+6u26QDr9hq68PLEjt/G9A+nWFbhLmT9vqNRmzZGk0aUBa/dAdDR3cHK2pUORRSaG+aWkZ7iGtBysP7Zbcdq0nDENRWxt6pft8+w4E9/JaN7KsfSnqKb2P7hFMsq3IXUH+84OYA20nQ9jSjra+2OWFWcm8FHp47iyXW13H7x6eRmhH7qVHtbGJKZSoF2t3XE6IIszps4jCfW1vK1iyY6vsCRz2f4ytLV7D9cQmPqn2lOXQ7dsHjjYr5//vcZmTvS0fjizcnJC6uPRqU0ryUNFbKb57s5fqKLZ9bX9ut+Hm8r5cU52jvGQQtmlXLwWDsrHF7Vz+czfOeZd3jhbS/NactoSlt68rZYr6aNVROG55GXkRq1dg1NGipk00sLmFqSz5LK/i0HW61TojvuwjNHMDQnnWVvOVdF5fMZvvfcZpa+VUNm/r9pSHlv43c8VNPGohSXML2sIGpJQ6unVMhEhJvnl3P7E5tYucfL2eP7XlGws9tH7dE2rpg2OgoRqt6kp7r4xMwxLH7TQ/3xE1Gfdt8Yww+Wb+ax1fv4jwtO45sXX4bI3VGNIZFVuAv5v3/u4lh7J0MivMiZljRUv1w+dRRFOeksCXH229qjbXT7jJY0YsCC2aV0+QxP97N6cbCMMfz3X7fy51X7+MJ54/jmxadrVWWYVbgLMQYu+tPnI94LTZOG6pfMtBSun13KK9sOUXu0794aHru7rY4Gd9744XlUuAtZ9lb0VvUzxrDo+W08tNLD588Zyx2XnqEJIwKseesMu+oi3y6kSUP1243z3AAhLQdbXa9jNGLJgtml7DnSEpX6b2MMP3lxGw++uZfPnF3Odz96piaMCDneWU+ny0O67/SIDxbWpKH6bUxBFh+eNIKla/b1OaeRx9tKTnoKxbmxM+V7MvvoB0aRk54S8QZxYww/+9t2Hnh9Lwvnu/nB5ZM0YUTQohWLaE37O20payPeC02ThhqQhfPLOdrayV83nXqRn2pvi3a3jSE5Gal8bNponn+7juYIrepnjOHul3dw32tVfGpeGT+8YrK+/xFU11zH4o2LaXT9lebU5XR0d0S0tKFJQw3I/NOGMmF4LksqPaesH6/2tuq64DFmwexS2jq7ef7tuogc/1f/2Mnv/r2HG+aU8aMrpmjCiLD3TE1ki2RpQ5OGGhAR4eazytm8/xgbahp73Ker20fN0VZdFzzGTC8tYOKI3IhUUf36lZ3c86/dLJhVyl1XTdHlW6Mg2lMT6TgNNWCfmDGG/3lpOw+v9JycojnQgcZ2OruNljRijLWqXxmLnt/KjoPNnD4yLyzH/c0/d/HrV3ZxTUUJP/3EBzRhREm0pybSkoYasJyMVK6uKOGFd+o43Nz+vts9J2e31ZJGrPn4jDGkpUjYShv3vrqbX/xjJ5+YMYafXz1VE0YC06ShBuXm+W46uw1L17z/y8c/JXp5sZY0Yk1RTjofmTSSZzbUcqJrcKv63ffaHu5+eQdXTR/N3ddOc3xCRBVZmjTUoIwblst5E4fx6OpqOrvf2xjn8baSmeZieF50p6xQoVkwu5SjrZ38Y+uhAR/jj69X8dOXtvOxaaP5X00YScGRpCEiRSLyDxHZZf99X4W4iEwXkUoR2SIib4vIAidiVX1bON/NoWMn+PuW9375VHtbKB+q3W1j1TnjixlTkDXgKqo/vbGXH7+wjY9+YBS/um4aqSn6GzQZOPUu3wH80xgzAfinfT1YK3CzMWYycAnwaxEpiGKMKkQXnD6c0qKs9y0H69HutjHN5RKunVXCG7vrQ5oSJtCSlR4WPb+VS6eM5NfXT9eEkUSceqevBJbY/y8BrgrewRiz0xizy/7/AHAYGBa1CFXIUlzCTfPcrNnbwLa6Y4C1Mts+bytuXRc8pl07qxSAJ9aGPonhI6uquXP5Fj4yaQT33DCDNE0YScWpd3uEMcY/suggMOJUO4vIHCAd2NPL7beKyFoRWXvkiLOLzCSr62aVkpnm4uFKDwAHj7XT0e3TkkaMG1OQxbkThvHkulq6fX1PYvjY6n18/9nNXHTmCH77yZmaMJJQxN5xEXlFRDb3cLkycD9jDSfu9WwVkVHAI8BnjAka9vjuMe43xswyxswaNkwLI04oyE7nquljeGbDfppaO/HUa3fbeLFgVin7G9t4Y3f9KfdbumYf33nmHT50xnDuvXEG6amaMJJRxAb3GWMu6u02ETkkIqOMMXV2Ujjcy35DgBeA7xpjVkUoVBUmN813s/StGp5YV0NWegqAljTiwEWThlOYncbjb9Vw/sSef3Q9vraGbz/zDudPHMbvbpxJRmpKlKNUscKpnwrLgYX2/wuB54J3EJF04BngYWPMk1GMTQ3Q5NH5zC4v5OHKavYeaSE91cXIIZlOh6X6kJGawidmlvD3rQfxHj/xvtufWlfLt556m3PGF3PfTRVkpmnCSGZOJY2fAR8WkV3ARfZ1RGSWiPzR3uc64Dzg0yKy0b5MdyZcFaqb55ezr6GVJ9fXYFyHOdw68DEAKnoWzC6ls9vwzIb979n+7Ib93P7kJs4+rZgHbp6lCUM5kzSMMV5jzIXGmAnGmIuMMQ329rXGmM/b///ZGJNmjJkecNnoRLwqdJdMGcnwvAwaW7s41r0n4quIqfCYOCKPGWUFPL723VX9ntu4n68/vpH544ZqwlAnaUuWCqu0FBdXzCgCoEP2R3wVMRU+C2aVsvPQceb/4Qb+vGYr/7lsI7PLi/jjwlkn26iU0qShwm7viUfw0UyHa0fEVxFT4XP5tNGkuLrYt28+3396D7PcRTz46dlkp+tk2OpdmjRUWNU11/HY1j9Qk3kjrSlvRnwVMRU+zR1HOJ7yGulmPCdSdsEyIQIAAAaHSURBVPLTa8vIydCEod5Lk4YKq5OriMm7Q2q0tBEfFq1YRHPa4zSlPklD5iJ+ueonToekYpAmDRVW0V5FTIWHf53pNrOfxrSHOOFr0hKi6pGWPVVYRXsVMRUep1pn+t6P3utQVCoWaUlDKaUlRBUyLWkopbSEqEKmJQ2llFIh06ShlFIqZJo0lFJKhUyThlJKqZBp0lBKKRUy8c9omShE5AhQHaHDFwOnXt4s/uhzig/6nGJfvD8ftzGmz6VPEy5pRJKIrDXGzHI6jnDS5xQf9DnFvkR7Pr3R6imllFIh06ShlFIqZJo0+ud+pwOIAH1O8UGfU+xLtOfTI23TUEopFTItaSillAqZJg2llFIh06TRTyKySETeFpGNIvJ3ERntdEyDJSJ3i8h2+3k9IyIFTsc0WCJyrYhsERGfiMRtN0gRuUREdojIbhG5w+l4BktEHhSRwyKy2elYwkVESkXkVRHZap9zX3U6pkjSpNF/dxtjphpjpgPPAz9wOqAw+AcwxRgzFdgJfNvheMJhM/AJYIXTgQyUiKQA9wKXApOAG0RkkrNRDdpDwCVOBxFmXcA3jDGTgHnAlxLgfeqVJo1+MsYcC7iaA8R9TwJjzN+NMV321VVAiZPxhIMxZpsxZofTcQzSHGC3MabKGNMBLAWudDimQTHGrAAanI4jnIwxdcaY9fb/zcA2YIyzUUWOLsI0ACJyF3Az0AR80OFwwu2zwDKng1CA9cVTE3C9FpjrUCwqBCJSDswAVjsbSeRo0uiBiLwCjOzhpu8aY54zxnwX+K6IfBv4MnBnVAMcgL6ek73Pd7GK2o9GM7aBCuU5KRUtIpILPAV8LahGIqFo0uiBMeaiEHd9FHiROEgafT0nEfk0cDlwoYmTwTv9eJ/i1X6gNOB6ib1NxRgRScNKGI8aY552Op5I0jaNfhKRCQFXrwS2OxVLuIjIJcA3gSuMMa1Ox6NOeguYICJjRSQduB5Y7nBMKoiICPAnYJsx5pdOxxNpOiK8n0TkKeB0wIc1Bfttxpi4/vUnIruBDMBrb1pljLnNwZAGTUQ+DvwGGAY0AhuNMRc7G1X/ichlwK+BFOBBY8xdDoc0KCLyF+ACrGnEDwF3GmP+5GhQgyQi5wCvA+9gfS8AfMcY86JzUUWOJg2llFIh0+oppZRSIdOkoZRSKmSaNJRSSoVMk4ZSSqmQadJQSikVMk0aKumISIGIfDHg+gUi8nw/j/HpSM9wLCI/FJHbI/kYSvWXJg2VjAqAL/a516l9Goj7afGV6i9NGioZ/Qw4zV4T5W57W66IPGmvK/KoPcoXEakQkddEZJ2IvCwio0TkGmAW8Kh9jCwR+YGIvCUim0Xkfv/9/UQkX0SqRcRlX88RkRoRSRORW+z7bhKRp0QkOzhgEfm3f10QESkWEY/9f4q9Hspb9nooX7C3jxKRFXZ8m0Xk3Ai9lirJaNJQyegOYI8xZrox5r/sbTOAr2GtWzEOONueT+g3wDXGmArgQeAuY8yTwFrgRvsYbcBvjTGzjTFTgCysebxOMsY0ARuB8+1NlwMvG2M6gaft+07Dmlb7c/14Lp8Dmowxs4HZwC0iMhb4pH386cA0+7GVGjSdsFApyxpjTC2AiGwEyrGmH5kC/MMuOKQAdb3c/4Mi8k0gGygCtgB/DdpnGbAAeBVrHqnf2duniMiPsarNcoGX+xH3R4CpdukHIB+YgDVv1YN24nvWGKNJQ4WFJg2lLCcC/u/G+mwIsMUYM/9UdxSRTKwEMMuY/9/e3YJEEIRhHP8/zXIIgtlmsWm6InbLBZtdMBiNgskkoibhLhw2tYkiFptVECw2waBymhQEP3gNM6Keezpg9PmVhd3ZYTbsPsy7MBOXkhaBvoqmu8CSpAFgDDjK59tAIyJO82rDExX3vvBRGfjct4C5iPgWNJLGgUmgLWklIjZ/eg6zEi5P2X90D9QK2p0Dg5LqkJa/ljRS0cf7R/w276kwRYWIeCDNANaAvYh4zZdqwFWeFUz3GMsFKWjo6v8QmM33Imk4/y8ZAm4iogm0gNGC5zX7lWca9u9ExJ2kY0lnwAGw36PdUy77rEvqJ70vq6TSUxvYkPQI1IEmaV/ya1Iw9LIF7PB1NrFA2umtk49VgbYMbEua6Rpvi1RKO8k/3ztAI/c/L+kZeCDtNGn2Z17l1szMirk8ZWZmxRwaZmZWzKFhZmbFHBpmZlbMoWFmZsUcGmZmVsyhYWZmxd4AjxOiF5DORrgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "diff = []\n",
    "for i in range(len(thetas)):\n",
    "    diff.append(values_physical[i] - values2_physical[i])\n",
    "\n",
    "plt.plot(thetas, diff, 'g^')\n",
    "plt.plot(thetas, diff)\n",
    "plt.title('difference between simultaneous and naive for Tokyo')\n",
    "plt.xlabel('theta values')\n",
    "plt.ylabel('energy sums')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the simultaneous measurement circuits\n",
    "more automations can be done in here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the number of shots\n",
    "shots = 1000\n",
    "\n",
    "def measure_zi_and_iz(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    zi = (counts.get('00', 0) + counts.get('10', 0) - counts.get('11', 0) - counts.get('01', 0))/shots\n",
    "    iz = (counts.get('00', 0) + counts.get('01', 0) - counts.get('10', 0) - counts.get('11', 0))/shots\n",
    "    return zi, iz\n",
    "\n",
    "def measure_xx_and_yy(theta):\n",
    "    circuit = get_ucc_ansatz(theta)\n",
    "    circuit.h(1)\n",
    "    circuit.cx(1, 0)\n",
    "    circuit.cz(0, 1)\n",
    "    circuit.h(0)\n",
    "    circuit.h(1)\n",
    "    circuit.measure(0, 0)\n",
    "    circuit.measure(1, 1)\n",
    "    job = qiskit.execute(circuit, backend=backend, shots=shots)\n",
    "    job_monitor(job)\n",
    "    counts = job.result().get_counts(circuit)\n",
    "    xx = (counts.get('00', 0) + counts.get('10', 0) - counts.get('11', 0) - counts.get('01', 0))/shots\n",
    "    yy = (counts.get('00', 0) + counts.get('01', 0) - counts.get('10', 0) - counts.get('11', 0))/shots\n",
    "    return xx, yy\n",
    "\n",
    "def measure_simultaneously_hamiltonian(theta):\n",
    "    xx, yy = measure_xx_and_yy(theta)\n",
    "    zi, iz = measure_zi_and_iz(theta)\n",
    "    return 5.9 + .22 * zi - 6.1 * iz - 2.14 * xx - 2.14 * yy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the experiment with different theta value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-3.141592653589793\n",
      "the current values[] array is: \n",
      "[11.78468]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-2.6179938779914944\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-2.0943951023931957\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-1.570796326794897\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-1.0471975511965983\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-0.5235987755982996\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "-8.881784197001252e-16\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "0.5235987755982978\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "1.0471975511965965\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "1.5707963267948948\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001, 2.16244]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "2.094395102393194\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001, 2.16244, 5.65808]\n",
      "Job Status: job has successfully run\n",
      "Job Status: job has successfully run\n",
      "theta is: \n",
      "2.617993877991493\n",
      "the current values[] array is: \n",
      "[11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001, 2.16244, 5.65808, 9.46308]\n"
     ]
    }
   ],
   "source": [
    "import qiskit\n",
    "import numpy as np\n",
    "values2 = []\n",
    "for theta in np.arange(-np.pi, np.pi, np.pi / 6):\n",
    "    values2.append(measure_simultaneously_hamiltonian(theta))\n",
    "    \n",
    "    # print out the values after each runs in order to save progress from program collapse/network issue\n",
    "    print('theta is: ')\n",
    "    print(theta)\n",
    "    print('the current values[] array is: ')\n",
    "    print(values2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot the results for the simultaneous measurement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNX5x/HPk52sLAlrgLAGAREkgIgKKK5117aIC26IVavWrdr+7CKttWJV3EER3FeqtRZlU0EEgQBhhyRsSSBAWBLCkkwyeX5/ZLBIQYYkM3eW5/163ReZm5t7vwO88sw599xzRFUxxhgTviKcDmCMMcZZVgiMMSbMWSEwxpgwZ4XAGGPCnBUCY4wJc1YIjDEmzFkhMMaYMGeFwBhjwpwVAmOMCXNRTgfwRmpqqmZkZDgdwxhjgsrixYt3qmra8Y4LikKQkZFBdna20zGMMSaoiMhmb46zriFjjAlzVgiMMSbMWSEwxpgwZ4XAGGPCnBUCY4wJc1YIGlBxeTGDJw9m275tTkcxxhivWSFoQGPmjGFuwVzGzB7jdBRjjPGaFYIGsm5HIe8t2kSUuwOTciZZq8AYEzSsENSDqrJ48x7u+zCHC59dSkrFnbSsfJq4ykt47BtrFRhjgkNQPFkcaPZWVPHp0i28u6CAtdvKiY+JYF/kLPZGzSSp+lKSXTfwyfwc7upXSPcWbZ2Oa4wxP8kKwQlYXlTKO98X8NmyrRysctOjdTKPX3Ey32wbyxsrJuByu6iMWE2F+0KaVo3iihcXMnlkPAM7NXM6ujHGHJMVguPYV1nNZzlbeXfhZlZu2Uuj6EguPaU1Iwa0o1d6CiLC2PFzcbldtT8gsC/qCyoj1pJe8weufe17fn12F+4+pwuREeLsmzHGmKOwQnAMq7aW8e6CAj5duoX9LjfdWibx2GU9uLxPG5Ljon907NLRS496jv2V1Tz6r5WMm5XH9xt2MW54H1qmxPkjvjHGeE1U1TcnFnkduBjYoao9PfvGApcALmA9cJOqlh7vXFlZWeqP2UcPuKr5fFkx7ywsYFlhKbFREVzcq/bT/6ntGiNSt0/0Hy8u4tFPV9IoJpJ//OIUhmY2b+Dkxhjzv0RksapmHfc4HxaCs4B9wJuHFYLzgK9UtVpE/g6gqr893rl8XQjWbSvn3QWb+efSLZRXVNO5eSIj+rfjqlPTSYmPPv4JvJC/Yx93vbuEtdvKGX1WRx44P5PoSBu0ZYzxHW8Lgc+6hlR1johkHLFv+mEvvweu9tX1j6eiys3UFcW8u6CA7M17iImM4MKTWzKifzv6d2ha50//x9K5eSKf3jmIMZ+vZvycDSzctJvnhvehbdP4Br2OMcacKJ+1CAA8heDzQy2CI773b+ADVX37GD97G3AbQLt27fpu3uzV+go/UlxezPApw/ng6g9omdgSqP1k/u6CAqYsKaLsYBUdUhNqP/33TadpQswJX6Mu/rO8mIenLEcEnry6Fxf0bOWX6xpjwovjXUOeEBkcpRCIyO+BLOBK9SJAXbuG7vjPHYxfPJ5Rfe7gwrYP8+6CAhZs3E10pHBej5Zc278dAzs1a/BP/94o2HWAu95bwvKiMm4Y2J7fXXQScdGRfs9hjAldAVsIRORGYDRwjqoe8OY8dSkExeXFdHn2dKIrh5LkPpcITaZt00Zc078dP+/blrSk2BM6ny+4qmt48su1vDZ3I91bJfPCiD50TEt0OpYxJkR4Wwj8erdSRC4AHgIu9bYI1NWYOWOIr7yM5OrLqYxYRf+eC5n9wFDuGNI5IIoAQExUBP93cXcmjsxia9lBLnl+Lp8u3eJ0LGNMmPHlqKH3gCFAKrAd+CPwCBAL7PIc9r2q3n68c51oi6C4vJiOz3WkypUIUoNbdtMoqhEb7tnww72CQLO19CD3vL+URZv28PO+6fz5sh7Ex9hjHsaYunO8RaCq16hqK1WNVtV0VZ2oqp1Vta2q9vZsxy0CdTFmzhhqtAZ3xE7cshsAt7oDenro1o0b8d6o0/j12Z35eEkRl77wHeu2lTsdyxgTBkJyIPv8ovn/nfLBw+V2Ma9onkOJvBMVGcH952Xy1s0DKD1QxaUvzOW9hQX48j6OMcb49GZxQ/HXk8WBpKS8kvs+zOHbvJ1cckprHr+iJ0lxDfNwmzEmPDjeNWTqJy0pljdu6s+D52cydUUxFz8/lxVFZU7HMsaEICsEASwiQrhzaGfev+00XNU1XPnyd7w+d6N1FRljGpQVgiDQL6MpU+8+k8Fd03js89WMenMxpQdcFJcXM3jyYFsW0xhTL1YIgkSThBhevSGLRy/uzuzcHVw07lt+8+8XmFswN6BHQxljAp8VgiAiItxyRgem/Op0RGqYv3wAsdVZTMqZZK0CY0ydWSEIQr3SG9P9pH9TFbGRVNeDSHVbaxUYY+rMCkEQKi4v5u2Vr7Ij5jFq2E+Tiod5Y8kn1iowxtSJFYIg9MOT07KbktgxRGgyyQcf4E9f/9XpaMaYIGSFIAgd/uS0K2I9O2P+QUxNJl/ltLKhpcaYE2azmgWhpaOX/s++F7/OZ+y02j/vOruLA6mMMcHKCkGIuGNIJ/K2l/PU9Fw6pSVy4cm26pkxxjvWNRQiRIQnrurFqe0a85sPc2w6CmOM16wQhJC46EjGX59Fs4RYbn1zEdv3VjgdyRgTBKwQhJi0pFheG5lFeUU1t76RzUGX2+lIxpgAZ4UgBJ3UKpnnhvdh5dYy7v8oh5oaG0lkjDk2KwQhalj3FjxyYTemrtjGs7PynI5jjAlgNmoohI06syN52/fx3Kw8OqUlcFnvNk5HMsYEIGsRhDAR4a9XnEz/Dk158OPlLC3Y43QkY0wA8lkhEJHXRWSHiKw8bF9TEZkhInmeP5v46vqmVkxUBK9c15cWybGMenMxW0sPOh3JGBNgfNkimAxccMS+h4FZqtoFmOV5bXysaUIMr4/sR2WVm1veyGZ/ZbXTkYwxAcRnhUBV5wC7j9h9GfCG5+s3gMt9dX3zY11aJPH8iD6s27aXez+wkUTGmP/y9z2CFqpa7Pl6G9DiWAeKyG0iki0i2SUlJf5JF+KGZDbn0Yu7M2P1dp6cts7pOMaYAOHYzWKtnSbzmB9LVXWCqmapalZaWpofk4W2G0/PYMSAdrwyez0fLy5yOo4xJgD4uxBsF5FWAJ4/d/j5+mFPRPjzpT04vVMzHvnnchZtOrL3zhgTbvxdCD4DRnq+Hgn8y8/XN0B0ZAQvXXsq6U3iGf3WYgp3H3A6kjHGQb4cPvoeMB/IFJEiEbkFeAI4V0TygGGe18YBjeNjmDgyi2p3Dbe8sYjyiiqnIxljHOLLUUPXqGorVY1W1XRVnaiqu1T1HFXtoqrDVNX6JRzUMS2Rl6/ry/qS/dz93lLcNpLImIBSXF7M4MmDfb4euT1ZHOYGdU7lz5f24Ot1JTw+dY3TcYwxh3ls9hjmFsxlzOwxPr2OFQLDdae158bTM5g4dyPvLSxwOo4xBlhcuInPvxtATHVPJuVM8mmrwAqBAeD/fnYSZ3VN49FPVzJv/U6n4xgT9u75+AsiNIWqiC241e3TVoEVAgNAVGQEL4zoQ0ZqAr96ewkbd+53OpIxYSu7YBOF29tQHvklbtmFy+3yaavACoH5QXJcNBNHZhEhcMsbiyg7aCOJjHHCvVO+BNzsjf7oh32+bBVYITA/0r5ZAq9c15fC3Qe4690lVLtrnI5kTFjZvGs/Rdtbe1oD/x1Y6XK7mFc0zyfXtIVpzP8Y0LEZf738ZB6aspzHPl/NY5f1dDqSMWHj+a/yiY2KZu1D42iePN4v17RCYI7qF/3akl+yjwlzNtC5eSI3DMxwOpIxIW/Tzv18snQLIwdm0Dw5zm/Xta4hc0y/vaAb53Rrzp//vZpPl63zy4MtxoSz577KIzpSuH1IR79e1wqBOabICGHcNX3o0jyRBz5cw/ebNvr8wRZjwtWGkn18unQL1w1oT/Mk/7UGwAqBOY7E2Cgev7odlTX7Sa38HZOWvmmtAmN84Pmv8omJimD04E5+v7YVAnNcry17krLYF4nWtsRWDbFWgTENbH3JPv6Vs4UbBmaQlhTr9+tbITA/qbi8mEk5kyhnPhURa0iovJpJS9+2VoExDei5WXnERkVy21n+vTdwiBUC85PGzBlDjdaAQGnUm0SRSpxrmLUKjGkg+TvK+WzZVm44vT2pif5vDYAVAnMc84vm43K7AKiMXMHBiBwSXFfyXUG2w8mMCQ3jZuXTKDqS0Wf5/97AIfYcgflJS0cv/fHrgj1c8dI8bs1826FExoSO3O3lfL58K7cP7kTThBjHcliLwJyQPu2aMOyk5oyfs4GyAzYXkTH1MW5WHvHRkdx2pjP3Bg6xQmBO2H3nZlJeUc2Eb9c7HcWYoLVuWzlTVxRz46AMmjjYGgArBKYOurdO5uJerZj03SZ27qt0Oo4xQWncrFwSYqIY5XBrABwqBCLyGxFZJSIrReQ9EfHvY3Sm3n5zblcqqty89LW1Cow5UWuK9zJ1xTZuGpRB43hnWwPgQCEQkTbA3UCWqvYEIoHh/s5h6qdTWiJXnZrO2ws2U1x20Ok4xgSVcTPzSIqN4tYznG8NgHNdQ1FAIxGJAuKBrQ7lMPVw9zldUFWem5XvdBRjgsaqrWV8uWobN53RgZT4aKfjAA4UAlXdAjwFFADFQJmqTvd3DlN/bZvGc03/dnyUXcjmXba0pTHeGDczj6S4KG45o4PTUX7gRNdQE+AyoAPQGkgQkeuOctxtIpItItklJSX+jmm8dNfQzkRFCs/OzHM6ijEBb+WWMqav3s4tZ3QgpVFgtAbAma6hYcBGVS1R1Srgn8DpRx6kqhNUNUtVs9LS0vwe0nineXIcIwdm8GnOFnK3lzsdx5iA9uzMPJLjorg5gFoDcIKFQEQiRCS5ntcsAE4TkXgREeAcYE09z2kcdPvgTiTERPH09FynoxgTsFYUlTFzzXZuPbMjyXGB0xoALwqBiLwrIskikgCsBFaLyIN1vaCqLgA+BpYAKzwZJtT1fMZ5TRJiuOWMDny5ahsrisqcjmNMQHp2Zi4pjaK5aVCG01H+hzctgu6quhe4HPiC2r796+tzUVX9o6p2U9Weqnq9qtpTSUHu1jM70Dg+mqemr3M6ijEBZ1lhKbPW7mDUmR1ICrDWAHhXCKJFJJraQvCZp19ffRvLBJukuGhuH9yJ2bklLNq02+k4xgSUZ2fm0jg+mpGnZzgd5ai8KQTjgU1AAjBHRNoDe30ZygSnkZ7VlcZOW4eqfVYwBmpn7P16XQmjzuwYkK0B8KIQqOpzqtpGVS/SWpuBoX7IZoJMo5hI7hramYUbd/Nt3k6n4xgTEJ6dmUeTAG4NgBfrEYhIY+AGIOOI4+/2USYTxIb3b8uEORt4avo6zuySSu3AMGPC0+LNe5idW8JvL+hGYmzgLv/iTdfQVGqLwApg8WGbMf8jNiqSe87pwvKi2gdnjAlnz87MpWlCDDcMbO90lJ/kTYmKU9X7fJ7EhIwrT23DK7PX8/T0XIad1ILICGsVmPCzeHNtF+kjF3YjIYBbA+Bdi+AtERklIq1EpOmhzefJTNCKiozg3nO7ss6zDJ8x4eiZGXmkJsZwfYC3BsC7QuACxgLz+W+3kK1cbn7SxSe3olvLJJ6ZkUuVu8bpOMb41aJNu5mbv5PRZ3UiPiawWwPgXSG4H+isqhmq2sGzBcYk2iZgRUQI95+XyaZdB5iyuMjpOMb41TMzcklNjOW60wK/NQDeFYJ84ICvg5jQM+yk5pzStjHPzcqjstrtdBxj/GLBhl3MW7+L2wd3pFFMpNNxvOJNIdgP5IjIeBF57tDm62Am+IkID56XydayCt5dUOB0HGP84pmZuaQlBU9rALwbNfSpZzPmhA3q3IzTOjblxa/z+WW/tkHRX2pMXc1fv4vvN+zmDxd3Jy46OFoD4EUhUNU3/BHEhCYR4cHzM7nq5flMnreJO4Z0djqSMT6hqjwzM5fmSbGMGNDO6TgnxJtpqDeKyIYjN3+EM6Ghb/umDM1MY/zsDZQdrHI6jjE+MX/9LhZu3M0dQzoFVWsAvLtHkAX082xnAs8Bb/sylAk995+XSdnBKiZ+a58hTOg51BpomRzH8P7B1RoA7yad23XYtkVVnwV+5odsJoT0bJPCRSe3ZOLcjezaZ8tPmNDyXf4uFm3awx1Dg681AN51DZ162JYlIrfj3U1mY37kvnO7crDKzSuz1zsdxZgGc6g10Coljl/2a+t0nDrx5hf6Pw77upratQl+4ZM0JqR1bp7E5X3a8Ob8zdx6ZkdaJMc5HcmYevs2byeLN+9hzOU9iY0KvtYAeNc1NPSw7VxVHaWqth6hqZN7z+mKu0Z5/qs8p6MYU2+HWgOtU+L4RVa603HqzJuuoXs8i9eLiLwmIktE5Dx/hDOhp12zeH7Zry0fLCqkcLc9sG6C2+zcEpYWlHLn2Z2DtjUA3o0autmzeP15QDNqF65/oj4XFZHGIvKxiKwVkTUiMrA+5zPB5ddndyFChHGzrFVggldtayCPNo0b8fO+wXlv4BBvCsGhyeQvAt5U1VWH7aurccCXqtoNOAVYU8/zmSDSMiWO609rzz+XFJG/Y5/TcYypk2/WlbCssJQ7h3YmJsqbX6WBy5v0i0VkOrWFYJqIJAF1nldYRFKAs4CJAKrqUtXSup7PBKdfDelEo+hInpmZ63QUY07YoXsDbRo34uq+wXtv4BBvCsEtwMNAP1U9AMQAN9Xjmh2AEmCSiCz13HdIqMf5TBBqlhjLzWd04D/Li1m1tczpOMackK/W7mB5URm/Pjv4WwPg3aihGlVdcuhTu+fBsuX1uGYUcCrwsqr2oXZ204ePPEhEbhORbBHJLikpqcflTKC69cyOJMdF8fR0axWY4KGqPDszj7ZNG3FVCLQGwLsWQUMrAopUdYHn9cfUFoYfUdUJqpqlqllpaWl+DWj8I6VRNKMHd2LW2h0sKdjjdBxjvDJzzQ5WbCnj10O7EB0Z/K0BcKAQqOo2oFBEMj27zgFW+zuHCQw3DcogNTGGp6bZoykm8G3du5U7P/gPbRrHcMWpbZyO02C8eY7gHyLSo4Gv+2vgHRFZDvQGHm/g85sgER8TxR1DOjNv/S7m5e90Oo4xP+muT17FVdmCZs0XhUxrALxrEawBJojIAhG53TPqp15UNcfT7dNLVS9XVesXCGMjBrSjVUocY6evQ1WdjmPMUW3Zu5UF65pRJVuZufVxtu3b5nSkBuPNzeLXVHUQcAOQASwXkXdFZKivw5nwEBcdyd3ndGFpQSlfrd3hdBxjjuqOjycRXdOB0qh3cFPFmNljnI7UYLxq24hIJNDNs+0ElgH3icj7PsxmwsjVfdNp3yyep6bnUlNjrQITWDbt2cKSvDZUSh4HIufgcruYlDMpZFoF3twjeAZYR+0DZY+ral9V/buqXgL08XVAEx6iIyP4zbCurCney9SVxU7HMeZHRn/wHpGaxp7o10FqP6i41R0yrQJvWgTLgVNUdbSqLjzie/19kMmEqUtOaU3XFok8PSOXanedH143pkHt2e8ityCDAxGLqIxc8cN+l9vFvKJ5DiZrON6sR7AMyBT50fRCZcBmVbVHQk2DiYwQ7js3k9vfXszk+at5K/9OPrj6A1omtnQ6mgljL3ydj2gjvr3nfjJb/snpOD7hTYvgJeB7YALwKjAf+AhYZ9NRm4Z2fo8W9EpP4R8z1jB38/ch0/Q2walg1wHenL+Jn/dtS2bLJKfj+Iw3hWAr0Mcz3LMvtfcFNgDnAk/6MpwJPyLCTWekcbAynvjqYSF1Q84En7HT19W2VM/r6nQUn/KmEHT1TD0NgKquBrqp6gbfxTLhbHrROCojVpFSNZyamihrFRhHLCss5d/LtjIqDJZV9aYQrBaRl0VksGd7ybMvFqjycT4TZorLi5m8bBK7oyYRRVNiKy+0VoHxO1Xl8alraJYQw+jBnZyO43PeFIKRQD5wr2fbANxIbRGwh8pMgxozZww1WoMrci0HIuaTUn0VNTWNrFVg/GrWmh0s2Libe4d1ITHWmzE1we0n36HnQbLXVPVa4B9HOcSWlzINan7RfFxuFwCl0W/SqvIFGlVeHjLD9Ezgq3bX8MSXa+mYmsDw/u2cjuMXP1kIVNUtIu1FJEZVXf4KZcLX0tFLf/T6wY+W8a9lV/P5L19wKJEJNx8trl1C9ZXr+obUxHI/xZt3uQH4TkQeFZH7Dm2+DmYMwL3n1o7WeHaGLV5jfG9/ZTVPz8glq30Tzu/Rwuk4fuNNIVgPfO45NumwzRifa9O4ETec1p4pS4rI217udBwT4l77diMl5ZU8ctFJHPEQbUg77l0QVf0zgIjEe9YsNsav7hzamQ8WFTJ22jom3JDldBwTonaUVzB+znou7NmSvu2bOB3Hr7yZdG6giKwG1npen+IZQmqMXzRJiGH04I5MX72dxZtt6QrjG+Nm5uGqruGhC7o5HcXvvOkaehY4H9gFoKrLgLN8GcqYI918RgdSE2P5+5drbfEa0+Dyd+zj/UWFXDugHR1SE5yO43de3RJX1cIjdrl9kMWYY4qPieKeczqzcONuvllX4nQcE2L+/uVaGnkWSApH3hSCQhE5HVARiRaRB6hdvtIYvxrevx3tm8Xz9y/X2uI1psEs3LibGau386shnWiWGOt0HEd4UwhuB+4E2gBbqF1s/k5fhjLmaKIjI7j/vEzWbivns2VbnY5jQsChqSRaJMdy86AOTsdxjDdrFu9U1WtVtYWqNlfV61R1V30vLCKRIrJURD6v77lM+Lj45Fb0aJ3MP2asw1Vti9eY+pm6Yhs5haXcf24mjWIinY7jGG9GDaWJyO9EZIKIvH5oa4Br34N1MZkTFBEhPHRBNwp3H+TdBZudjmOCmKu6hienrSWzRRJX9U13Oo6jvOka+heQAswE/nPYVmcikg78DHitPucx4emsLqkM7NiM57/KZ19ltdNxTJB6Z8FmNu86wMMXdSMyInweHjsabwpBvKr+VlU/VNUph7Z6XvdZ4CHA2vbmhIkIv72wG7v2u5j47Uan45ggtLeiiudm5TGoczOGdE1zOo7jvCkEn4vIRQ11QRG5GNihqouPc9xtIpItItklJTZc0PxY77aNubBnSybMWc+ufZVOxzFB5uVv1rPnQBWPXBheU0kcizeF4B5qi0GFiOwVkXIR2VuPaw4CLhWRTcD7wNki8vaRB6nqBM/ymFlpaVaxzf+6/7xMDla5eeHrfKejmCCytfQgr8/dyOW9W9OzTYrTcQKCN6OGklQ1QlXjVDXZ8zq5rhdU1UdUNV1VM4DhwFeqel1dz2fCV+fmifwiqy3vfF9A4W6bBst45+kZuajCA+dnOh0lYHgzakhE5DoRedTzuq2I9Pd9NGOO795hXRGBZ2baNNXm+FZv3cuUJUXcOCiD9CbxTscJGN50Db0EDARGeF7vA15siIur6jeqenFDnMuEp5Ypcdw4KINPlm5h7bb69FiacPDEl2tJjovmziGdnY4SULwpBANU9U6gAkBV9wAxPk1lzAn41eBOJMVGMfbLdU5HMQHs27wS5uSW8OuzO5MSH+10nIDiTSGo8qxdrFD7gBk27NMEkMbxMdw+pBOz1u5g0abdTscxAaimRnl86lrSmzTi+oHtnY4TcLwpBM8BnwDNReSvwFzgcZ+mMuYE3XR6B1okx/L3L2yaavO/Ps3ZwprivTx4fiaxUeE7lcSxeDNq6B1qH/76G1AMXK6qH/k6mDEnolFMJPec05XszXuYtWaH03FMAKmocvPUtHWc3CaFS3q1djpOQPJ2PYK1qvqiqr6gqjY/kAlIP89Kp0NqAk9OW4vbpqk2HpPnbWJrWQWPXNSNiDCfSuJYvCoExgSD6MgIHjgvk9zt+/hk6Ran45gAsGe/ixe/zufsbs05vVOq03EClhUCE1IuOrklvdJTeGZGLpXVtpBeuHv+q3z2V1bz8IXhtw7xibBCYEKKiPDbC7qxpfQgb39f4HQc46DNu/bz1veb+EVWW7q2SHI6TkCzQmBCzqDOqZzROZUXv86nvKLK6TjGIWOnrSMqIoLfnNvV6SgBzwqBCUm/vaAbu/e7eHXOBqejGAfkFJby+fJiRp3ZgRbJcU7HCXhWCExIOjk9hZ/1asVrczdSUm7TVIeTQ+sQN0uI4bbBnZyOExSsEJiQ9cB5mVRW1/DCV3lORzF+NGvNDhZu3M29w7qQGBvldJygYIXAhKwOqQn8sl9b3l1YQMEum6Y6HFS7a/jbF2vomJrA8P7tnI4TNKwQmJB2zzldiIwQ/jHDJqQLBx9mF7G+ZD8PXdCN6Ej79eYt+5syIa1Fchw3D+rAv3K2smprmdNxjA/tr6zmmZm5ZLVvwvk9WjgdJ6hYITAhb/TgTqQ0imbsNGsVhLJXv91ASXklj1xk6xCfKCsEJuSlNIrmjiGd+GZdCd9v2OV0HOMDO8ormDBnAxed3JK+7Zs4HSfoWCEwYWHk6Rm0SonjCZumOiQ9OzMPV3UND55vU0nUhRUCExbioiO5d1gXcgpLmb56u9NxTAMpLi/m9PFX8MHCAq4d0I4OqQlORwpKVghM2Ljq1HQ6pSUwdto6qt22yF4oGDNnDPkFPZGIKu4+p4vTcYKW3wuBiLQVka9FZLWIrBKRe/ydwYSnqMgIHjw/k/wd+/jnEpumOtgVlxfzdvb3xLtPozTyI6rY43SkoOVEi6AauF9VuwOnAXeKSHcHcpgwdH6PlvRu25hnZuZSUWXTVAez/5v5JCkVd1MtOyiP/owxs8c4HSlo+b0QqGqxqi7xfF0OrAHa+DuHCU+HpqkuLqvgrfmbnY5j6mjj7i18uag9oknsiPkLrpp9TMqZxLZ925yOFpQcvUcgIhlAH2CBkzlMeBnYqRmDu6bx4jf57LVpqoOOqjLi9alE1XRiZ8xTVEXUzjDrVre1CurIsUIgIonAFOBeVd17lO/fJiLZIpJdUlLi/4AmpD14fialB6oYP3u901HMCXp2Zh7FO1tTGvUGByO//2G/y+1iXtE8B5MFL0em5hORaGqLwDuq+s+jHaOqE4AJAFlZWTbw2zSonm1SuPSU1kycu5GRAzNobnOtzhRNAAAQDElEQVTWB4XPlm1l3Kw8ft43nSev/sieIG4gTowaEmAisEZVn/b39Y055P7zulLtVp6zaaqDwpKCPTzw0TL6ZzTlr1ecbEWgATnRNTQIuB44W0RyPNtFDuQwYa59swRGDGjH+wsLWbh5I4MnD7abjQFqS+lBbntzMS2T43jl+r7ERNkjUA3JiVFDc1VVVLWXqvb2bFP9ncMYgLvO7kx0ZAT3fDSTuQVz7WZjANpXWc0tkxdRWe3m9RuzaJoQ43SkkGNl1YS15klx/LJ/GsU7WxNd3cOGIAYYd41y7/tLyduxjxdHnErn5klORwpJVghM2NtaM5lqKSLN9TC4U61VEECe+GINM9fs4E+XdOesrmlOxwlZVghMWCsuL+atla+xPeYxIILGBx9m8tL3rVUQAD5YVMCr325k5MD2XD8ww+k4Ic0KgQlrY+aMoUZrqI7Yys6YJ4jWdJIO3s1j31irwEnz1+/i95+s5KyuaTx6sc1A42tWCExYm180H5fbBUBF5DJ2R48nzp3FrBXWF+2UjTv3c/vbi+mQmsALI/oQZWsP+5wjD5QZEyiWjl76P/se/XQlb30PH2UX8vOstg6kCl9lB6q4ZfIiIiOEiSP7kRwX7XSksGCl1pgj/OGS7gzq3IzffbKC7E27nY4TNqrcNdzx7mIK9xzglev60q5ZvNORwoYVAmOOEB0ZwUsj+pLeJJ7Rby2mcPcBpyOFPFXlj5+t4rv8Xfztyl7079DU6UhhxQqBMUeREh/Nqzdk4XLXMOrNbPZVVjsdKaRN+m4T7y4o4FdDOnF133Sn44QdKwTGHEPn5om8OOJU8nbs4973c6ipsbkPfeHrtTv4y39Wc36PFjx4XqbTccKSFQJjfsJZXdN49GcnMXPNdsZOX+d0nJCzbls5v35vKSe1SuaZX/YmIsImknOCjRoy5jhGnp5B7o59vPzNero0T+TKU63roiHs3FfJzZMXER8TycSR/YiPsV9HTrEWgTHHISL8+dIenNaxKQ9PWcHizbZIen1VVLkZ/dZidu2v5LWRWbRMsfUgnGSFwBgvREdG8PK1fWnVOI7Rb2WzpfSg05GClqry8JTlLN68h6d/0Zte6Y2djhT2rBAY46UmCTFMHJlFZVUNo97I5oDLRhLVxYtf5/NpzlYeOK8rF53cyuk4BisExpyQzs2TeH5EH9Zu28tvPrCRRCdq6opinpqeyxV92nDn0M5OxzEeVgiMOUFDMpvz+591Z9qq7Tw9I9fpOEFjeVEp932YQ9/2TfjblbbUZCCx2/TG1MHNgzLI3VbOC1/n06VFIpf1buN0pIBWXHaQW9/IJjUxlvHX9yUuOtLpSOYw1iIwpg5EhDGX96R/h6Y8+PFycgpLnY4UsA64qrn1jWwOuNxMHNmP1MRYpyOZI1ghMKaOYqIieOW6vrRIjmXUm9kUl9lIoiPV1Cj3vp/DmuK9PH9NHzJb2vTegciRQiAiF4jIOhHJF5GHnchgTENomhDDxJH9OFBZzag3bSTRkcZOX8f01dv5v591Z2i35k7HMcfg90IgIpHAi8CFQHfgGhGxJYhM0OraIonnrunDqq17eeCjZTaSyOPjxUW8/M16Rgxox02DMpyOY36CEy2C/kC+qm5QVRfwPnCZAzmMaTDnnNSCRy7sxtQV2xg3K8/pOI4qLi9mwEvX8vCU5Qzq3Iw/X9rDRggFOCcKQRug8LDXRZ59PyIit4lItohkl5SU+C2cMXU16syO/LxvOuNm5fH58q1Ox3HMI9OeYmvBxcTGlvPSiL5E21KTAS9g/4VUdYKqZqlqVlpamtNxjDkuEeEvV/SkX0YT7v9wGcuLwm8k0YJNG5m5pAsgFEQ8wsGaXU5HMl5wohBsAQ5fCDbds8+YoBcbFcnL1/UlNbF2JNH2vRVOR/KLbWUVPPLPFQx/ZSURNamUxDyOS7YwZvYYp6MZLzhRCBYBXUSkg4jEAMOBzxzIYYxPpCbG8trILMorqrntzWwqqtxOR/KZ0gMu/vbFGgaP/ZqPsgvZH/0FW+JupTJyBS63i0k5k9i2b5vTMc1x+L0QqGo1cBcwDVgDfKiqq/ydwxhfOqlVMuOG92H5ljIe/Hg5qqE1kuigy81L3+Rz1pNfM2HOBn52civO6DuLstiJ1EjZD8e51W2tgiDgyBQTqjoVmOrEtY3xl3O7t+Ch87vx9y/X0qV5Inef08XpSPVW5a7hw+xCxs3MY0d5JcNOas4D52fSrWUyfcbfhMvt+tHxLreLeUXzHEprvGVzDRnjQ7cP7kje9nKenpFLl+aJXBik0y7X1Cj/WVHMP6avY9OuA2S1b8KL155Kv4ymPxyzdPRSBxOa+rBCYIwPiQiPX3kyG3ft574Pl9G2aTw926Q4Hctrqsq3eTt5ctpaVm7ZS2aLJCaOzOLsbs3t2YAQErDDR40JFXHRkUy4Posm8dGMejObFVsLGDx5cMDfRM0pLGXEqwu44fWFlB6o4ulfnMLUe87knJNaWBEIMVYIjPGDtKRYXh2ZRemBKq6fNIe5mxcE7E3U/B37uP2txVz+4nfkbi/nT5d0Z9b9g7ny1HQiI6wAhCLrGjLGT3q0TuEPl7bnkSluUiN+z3uL5nFd940MaJ9BRAD8gt1aepBxM/P4aHEhjaIj+c2wrtxyZgcSY+3XRKizf2Fj/OjbHc9TFrONJNc1NKo4lWvGryYxNpeT26TQq20KvdMbc0rbxrRKifNb98ue/S5enr2eyfM2gcKNp3fgzqGdaGbrBoQNKwTG+ElxeTGTciZREVlBady/iNY2JEoPrurxZ/K3u3h97kaq3LXPG6QmxtK7bQq9PIXhlPQUGsfHNGieA65qXp+7kfGzN7DPVc2VfdK5d1gX2jaNb9DrmMBnhcAYPxkzZww1WlP7QmqokkL2R26nIr4j/7rrRSqr3awpLmd5USk5haUsKyxl5podP/x8+2bxnJLemF7pKfRu25gerVNoFHPiSz5WuWt4f2EB42bls3NfJcNOasGD52faojFhzAqBMX4yv2j+Tz5wFRsVSe+2jendtjE3DKz9/t6KKlYWlZFTVMrywjIWbdrNZ8tqZzaNjBC6tkj6b8shvTFdWyQSdcRsn8XlxQyfMpz3rnqfRetreHpGLpt3HaB/RlPGX38qfds3xYQ3CYZH37OysjQ7O9vpGMYEhB17K1hWVMaywlKWFdW2HPZW1K6MFhcdQc/WKZzS9r8thyfmP8ibixbRKeo37N2fQreWSfz2gm4MyUyzYaAhTkQWq2rWcY+zQmBMcFNVNu068KMupVVb91JZXdsNVcNBImiEW7Yz5tIB3DCge0CMUjK+520hsK4hY4KciNAhNYEOqQlc1rt2jacqdw3rtpXz8BcvsWBTMRWyHlfMNyzafSM3RrzobGATcOyBMmNCUHRkBM2SDzBr+2PsjH6efVFTcdUcsGmhzVFZITAmRP1olJKHTQttjsYKgTEh6nijlIw5xO4RGBOibFpo4y1rERhjTJizQmCMMWHOCoExxoQ5KwTGGBPmrBAYY0yYC4opJkSkBNjso9OnAjt9dG4nhNr7AXtPwcLeU+Bpr6ppxzsoKAqBL4lItjdzcQSLUHs/YO8pWNh7Cl7WNWSMMWHOCoExxoQ5KwQwwekADSzU3g/YewoW9p6CVNjfIzDGmHBnLQJjjAlzYV8IRGSMiCwXkRwRmS4irZ3OVF8iMlZE1nre1yci0tjpTPUlIj8XkVUiUiMiQT2KQ0QuEJF1IpIvIg87nae+ROR1EdkhIiudztIQRKStiHwtIqs9/+fucTqTr4V9IQDGqmovVe0NfA78welADWAG0FNVewG5wCMO52kIK4ErgTlOB6kPEYkEXgQuBLoD14hId2dT1dtk4AKnQzSgauB+Ve0OnAbcGQL/Rj8p7AuBqu497GUCEPQ3TVR1uqpWe15+D6Q7machqOoaVV3ndI4G0B/IV9UNquoC3gcuczhTvajqHGC30zkaiqoWq+oSz9flwBqgjbOpfMvWIwBE5K/ADUAZMNThOA3tZuADp0OYH7QBCg97XQQMcCiLOQ4RyQD6AAucTeJbYVEIRGQm0PIo3/q9qv5LVX8P/F5EHgHuAv7o14B1cLz35Dnm99Q2c9/xZ7a68uY9GeMvIpIITAHuPaLnIOSERSFQ1WFeHvoOMJUgKATHe08iciNwMXCOBskY4RP4dwpmW4C2h71O9+wzAUREoqktAu+o6j+dzuNrYX+PQES6HPbyMmCtU1kaiohcADwEXKqqB5zOY35kEdBFRDqISAwwHPjM4UzmMCIiwERgjao+7XQefwj7B8pEZAqQCdRQO8Pp7aoa1J/QRCQfiAV2eXZ9r6q3Oxip3kTkCuB5IA0oBXJU9XxnU9WNiFwEPAtEAq+r6l8djlQvIvIeMITamTq3A39U1YmOhqoHETkD+BZYQe3vBYDfqepU51L5VtgXAmOMCXdh3zVkjDHhzgqBMcaEOSsExhgT5qwQGGNMmLNCYIwxYc4KgQl6ItJYRO447PUQEfn8BM9xo69nnhWRP4nIA768hjF1YYXAhILGwB3HPeqn3QgE/RTkxtSFFQITCp4AOnnWlBjr2ZcoIh971mV4x/O0KCLSV0Rmi8hiEZkmIq1E5GogC3jHc45GIvIHEVkkIitFZMKhnz9ERFJEZLOIRHheJ4hIoYhEi8goz88uE5EpIhJ/ZGAR+ebQugoikioimzxfR3rWk1jkWU9itGd/KxGZ48m3UkTO9NHfpQlDVghMKHgYWK+qvVX1Qc++PsC91M753xEY5Jk/5nngalXtC7wO/FVVPwaygWs95zgIvKCq/VS1J9CI2nmbfqCqZUAOMNiz62JgmqpWAf/0/Owp1E5hfMsJvJdbgDJV7Qf0A0aJSAdghOf8vYFTPNc2pkGExaRzJiwtVNUiABHJATKonZqiJzDD8wE/Eig+xs8PFZGHgHigKbAK+PcRx3wA/BL4mto5g17y7O8pIn+htssqEZh2ArnPA3p5WikAKUAXaucoet1TzD5VVSsEpsFYITChqvKwr93U/l8XYJWqDvypHxSROGp/qWepaqGI/AmIO8qhnwGPi0hToC/wlWf/ZOByVV3mmQV2yFF+tpr/tsgPP7cAv1bV/ykeInIW8DNgsog8rapv/tT7MMZb1jVkQkE5kOTFceuANBEZCLVTDYtIj6Oc49Av5p2eOemv5ihUdR+1n9THAZ+rqtvzrSSg2PPp/dpjZNlEbfHgiPNPA37l+VlEpKvn/kN7YLuqvgq8Bpzqxfs1xivWIjBBT1V3ich3nsXTvwD+c4zjXJ4ul+dEJIXa///PUtvtMxl4RUQOAgOBV6ldJ3kbtb/sj+UD4CN+/Kn/UWpXtCrx/Hm0IvUU8KGI3HZE3teo7cZa4rlBXQJc7jn/gyJSBeyjdkU9YxqEzT5qjDFhzrqGjDEmzFkhMMaYMGeFwBhjwpwVAmOMCXNWCIwxJsxZITDGmDBnhcAYY8KcFQJjjAlz/w/ytO3QUljXCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "values2_physical = [11.78468, 12.558240000000001, 11.13888, 8.84964, 5.1592400000000005, 1.5702000000000012, -0.7365999999999993, -1.1382, -0.11680000000000001, 2.16244, 5.65808, 9.46308]\n",
    "thetas = []\n",
    "for theta in np.arange(-np.pi, np.pi, np.pi/6):\n",
    "    # put the theta values to a list for plotting\n",
    "    thetas.append(theta)\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(thetas, values2_physical, 'g^')\n",
    "plt.plot(thetas, values2_physical)\n",
    "plt.xlabel('theta values')\n",
    "plt.ylabel('energy sums')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot the difference between results obtained from the simulator and from Tokyo 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulated values copied from the other ipython notebook\n",
    "values2_simulated = [12.2842, 13.644680000000001, 12.96908, 10.20912, 6.76572, 2.5197200000000013, -0.45851999999999904, -1.7290799999999997, -1.2657199999999995, 1.5503200000000001, 5.36568, 9.352079999999997]\n",
    "diff_phys_sim = []\n",
    "for theta in range(len(thetas)):\n",
    "    diff_phys_sim.append(values2_physical[theta] - values2_simulated[theta])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX6x/HPkwoJgUAINUDoXVpEARVQVLChuPaGqyIi9rIqtiUWdm1YwIYGXRFhXQuKCiJdaiBBSqihJCEQCJCeTJI5vz8y+otIYIbM5M5MnvfrNa9k7ty59zsa8uTcc885YoxBKaWUclaA1QGUUkr5Fi0cSimlXKKFQymllEu0cCillHKJFg6llFIu0cKhlFLKJVo4lFJKuUQLh1JKKZdo4VBKKeWSIKsDeELjxo1NbGys1TGUUspnrFu37rAxJtqZff2ycMTGxpKYmGh1DKWU8hkistfZffVSlVJKKZdo4VBKKeUSLRxKKaVcooVDKaWUS7RwKKWUcokWDqWUZTLzMhk8fTAH8g9YHUW5QAuHUsoy8UvjWb5vOfFL4q2OolyghUMpZYnMvEw+XfczoWVxJCQnaKvDh2jhUEpZ4ukF/6Jh0bM0sT1LUGlfbXX4EC0cSqkal5G7nx8SmyAmHJvsJbL4Pj5d/6O2OnyEFg6lVI27fcYXhJb34mjwB2SF/BODnfpFj/D8whetjqacoIVDKVWjVqdms31vRwoCl5AfOI/ygCwOh7xKkL01C39rijHG6ojqFLRwKKVqTHZ+Cfd/kUTbxhGkPvMi5nmDec5QNHEdDw3rjC2/DzNW77M6pjoFLRxKqRphtxsemr2Bo4WlvHNjHyLqBP/p9fvP78jgTtFM/G4LyWnHLEqpnGFp4RCRj0UkS0Q2VfG6iMhbIrJTRH4Tkb41nVEp5R7vLd3F0u2HeOaybnRv0eAvrwcECJOv6010RCjjPlvHkQKbBSmVM6xucUwHhp/k9RFAR8djDPBuDWRSSrnZ2j1HeG3+di49ozk3n9W6yv0ahofw3s39OFxg44Evkii3a3+HN7K0cBhjlgJHTrLLSOBTU2EVECkizWsmnVLKHY4U2Lh/ZhIxDesyaVRPROSk+/eMacDEK7qzbMdhJi/YXkMplSusbnGcSksgrdLzdMc2pZQPsNsNj8xOJjvfxpQb+/6lX6Mq1/dvzbVxMby9cCe/pBz0cErlKm8vHE4TkTEikigiiYcOHbI6jlIK+HBZKou2HeLpy7rSo+Vf+zVOZuLIHnRvUZ+HZiWzL7vQQwnV6fD2wpEBtKr0PMax7S+MMR8YY+KMMXHR0U6tt66U8qB1e4/w73nbuKRnM245u43L768THMh7N/dDRBj72TqKS8s9kFKdDm8vHHOAWx13V50N5BhjMq0OpZQ6uaMFNu77PIkWkXWYdPUZp+zXqEqrRmFMvq43KQdyefqbTTo40EtYfTvuTGAl0FlE0kXkDhEZKyJjHbv8AKQCO4EPgXEWRVVKOckYw2NfbuBQfglTbuxLfSf7NaoytEsT7ju/I1+uS2fmmrRTv0F5XJCVJzfG3HCK1w1wbw3FUUq5wUfLd7MgJYvnLu/GGTGRbjnmAxd0JGnfUZ6fs5keLeu77bjq9Hj7pSqllA9J2neUST9u5eLuTRk9MNZtxw0MEN66vg/REaHc89l6jurgQEtp4VBKuUVOYSnjP0+iWYM6/PvqXqfdr1GVhuEhTL2pL4fySnhgVrIODrSQFg6lVLUZY3j0yw1k5RXzzo19aRBWvX6NqvRqFcnzV3Rn6fZDvPnLDo+cQ52aFg6lVLUl/LqHn7cc5B/Du9C7lWf7H27o34q/9YvhrV92sGhrlkfPpU5MC4dSqlo2pB3j5R9TGNa1KXec09bj5xMRXriyB92a1+fBWcmkHdHBgTVNC4dS6rTlFJVy7+fraRJRh1evOf3xGq76fXCgMUYHB1pAC4dS6rQYY/jHl79xIKeYt27oQ2RYSI2ev3VUGG9c15vN+3N59tsTrsygPEQLh1LqtHyyYg8/bT7A48M7069NQ0syXNC1KeOHdmB2Yjqz1urKgTVFC4dSymUb03N46YetnN+lCXee087SLA9d2IlzOzbmmW83szE9x9IstYUWDqWUS3KLK/o1ouqF8No1vQgIqJl+jaoEBghvXt+HxuEh3DNjHccKdXCgp2nhUEo5zRjDk//bSMaxIt65sQ8Nw2u2X6MqjcJDmHpzP7JyS3hwVjJ2HRzoUVo4lFJO+2zVXuZuzOTRizrTr00jq+P8Se9WkTx7eTcWbzvEWwt1cKAnaeFQSjllU0YO8d+nMKRzNHefZ22/RlVuOqs1o/q25M1fdrB4mw4O9BQtHEqpU8orLmX85+tpFB7C69f2trxfoyoiwotX9qRz0wge+EIHB3qKFg6l1EkZY3jyq42kHS3irRv60MhL+jWqUjckkPdv6YfdGMbNWK+DAz1AC4dS6qQ+X7OP73/L5OELO9G/rXf1a1SlTVQ4r1/bm40ZOTw/Z7PVcfyOFg6lVJW27M/ln99t4dyOjblncHur47jkwm5NGTekPV+sTWP2Wl050J20cCilTii/pIzxn68nsm4wb1znvf0aJ/PIRZ0Z1CGKZ77dxKYMHRzoLlo4lFJ/YYxhwtcb2ZNdwFs39KFxvVCrI52W31cObOQYHLgtK43B0wdzIP+A1dF8mhYOpdQfMvMyGTx9MB8s38y3yft5aFgnzm4XZXWsaomqF8rUm/pyIKeYmxMWsHzvr8Qvibc6lk/TwqGU+kP80nhW70njXz+mck6Hxowb2sHqSG7Rp3VDHhgWw6GjTYi03U1C0n+01VENWjiUUkBFayMh6XOiSh6n1OTy5GUtCPTBfo2qbC54j/zgr4kov4TIwud5av4rVkfyWVo4lFJARWujrm04waYVx0InM3XdJKsjuU1mXibTNySQHfQRh0JeJsjeigVr4/h2w3aro/kkLRxKKUdrYwbhtssoClhPvqwjITnBby7nxC+Nx27sABQG/kpm6EOUy1EemLmdyQu2U66TIrpEC4dSivil8dSxnU8gDckJmgVAuSn3m07klekrsZX//3TrZQEZZIY8TEi9ZCYv2MHohDUcKdDp2J0VZOXJRWQ48CYQCEwzxkw67vXRwCtAhmPTO8aYaTUaUqlaYEXaGsJLx1McsJmSwIqR1rZyGyvSV1iczD2S7k464XZjDDPXpPH8nM1c+tYyptzUl76trVnN0JeIMdY00UQkENgOXAikA2uBG4wxWyrtMxqIM8aMd+XYcXFxJjEx0Y1plfJvs9em8fj/fmP67WcypHMTq+PUuI3pOYz7fB0HcoqZcElXbhsYi4j/3BjgDBFZZ4yJc2ZfKy9V9Qd2GmNSjTE24AtgpIV5lKqVyu2Gd5fsomfLBgzuFG11HEv0jGnA9+PPZXCnaJ7/bgvjZyaRX1JmdSyvZWXhaAlUnkAm3bHteFeLyG8i8qWItKrqYCIyRkQSRSTx0KFD7s6qlN+auzGT3YcLuHdo+1r3V3ZlDcKC+eCWOJ4Y0YUfN2ZyxTvL2XYgz+pYXsnbO8e/A2KNMWcAPwOfVLWjMeYDY0ycMSYuOrp2/tWklKvsdsPURTvp0KQeF3VrZnUcywUECGMHt+fzu84mt6iMkVOW89X6dKtjeR0rC0cGULkFEcP/d4IDYIzJNsaUOJ5OA/rVUDalaoVftmax9UAe44a098lJDD3l7HZR/HD/OfSKieTh2Rt46uuNuq5HJVYWjrVARxFpKyIhwPXAnMo7iEjzSk+vAFJqMJ9Sfs0YwzuLdtKqUV2u6NXC6jhep0n9Osy486yKFsjqffztvRW6oqCDZYXDGFMGjAfmUVEQZhtjNovIRBG5wrHb/SKyWUQ2APcDo61Jq5T/+XVnNhvSjjF2cHuCAr39qrU1ggIDeGJEFz68NY692YVc+tYyFmw5aHUsy1l2O64n6e24Sp3aDR+sIvVwPksfH0poUKDVcbzevuxCxn2+jk0ZudwzpD2PXNjJrwqur9yOq5SyyLq9R1mZms1d57bTouGk1lFhfDl2IDf0b827i3dx80erycortjqWJbRwKFULTVm0k4Zhwdx4Vmuro/iUOsGBvDyqJ69d04vktGNc+tZyVqdmWx2rxmnhUKqW2bw/h4Vbs/j7oLaEhVg665DPurpfDN/cO4iI0CBunLaa95bswh8v+1dFC4dStczURbuICA3i1oGxVkfxaV2a1efb8YMY3r0Zk37cyl2friOnqNTqWDVCC4dStcjOrHx+2JTJLQPa0KBusNVxfF5EnWDeubEPz13ejcXbsrjs7WVsysixOpbHaeFQqhZ5d/EuQoMCuOOctlZH8Rsiwu2D2jLr7gGUlRtGvbuCmWv2sT93P4OnD/abNU0q08KhVC2RdqSQb5IzuKF/a6LqhVodx+/0a9OQ7+87h7PaNuLJrzZyzbQ5/Lp3jd+saVKZFg6laokPlqYSIDDmvHZWR/FbUfVCmX57f+44txnpWS1pUvwqMxJ/Y39eptXR3EoLh1K1QFZuMbMS0/hbvxiaN6hrdRy/FhggZPIRR+rEE0AYkUVPMOy1ZXy8fDd5xf7Rea734ilVC0xbvpuycjt3n9fe6ih+LzMvk4TkBIqlmPzQ9YSVDyDSdhUTvw/n9Z+3c01cDKMHxtImKtzqqKdNC4dSfu5ogY3PVu3l8l4tiG3su7+sfEX80njsxl7xRMopDFpOWega/tbhEaK5hv+s3Mv0FXu4oEsT/j6oLQPaR/ncOihaOJTycwkr9lBoK2fckA5WR6kVVqavxFZu+9M2W7mNLTk/knT3Szx5SVc+W7WXGav3sSBlNV2aRTB6YCxX9mlJnWDfmP5FJzlUyo/lFZcyaNJCzm4XxQe3OjV/naohxaXlzNmwn4+X72brgbw/poC55exYmjWoU+N5XJnkUFscSvmxGav3kVtcxvjztbXhbeoEB3JtXCuu6RfDqtQjJPy6m6mLd/H+klQu6dmc2wfF0qd1Q6tjnpAWDqX8VHFpOdOW7ebcjo05IybS6jiqCiLCgPZRDGgfxb7sQj5ZuYfZa9OYs2E/fVpHcvugtozo0YxgL5rC3XuSKKXcatbaNA7nlzB+qLY2fEXrqDCeuawbK5+6gOcv78bRAhv3z0zi3H8tYsqinRwtsJ36IDVA+ziU8kO2MjtDXllEi8i6/HfsAJ+7a0dVsNsNi7dn8fHyPSzfeZjQoACu6tOS2we1pXOzCLeeSxdyUqqW+yYpg/05xdx7fgctGj4sIEA4v0tTPrvzLOY/dB6j+sbwdVIGF09eyk3TVrFgy0Hs9oo//jPzMmtsbixtcSjlZ8rthmGvLyE8NJDvxp+jhcPPHC2wMXPtPv6zci+ZOcW0iQpj9MBYVh5+nY82TGVsv7FMuXSKy8f1WItDRAJEpL7LiZRSNWbuxkx2Hy7g3iHa2vBHDcNDGDekA0sfH8o7N/YhKjyEf363hZ9WnkeDkjtISPrM462OUxYOEflcROqLSDiwCdgiIo95NJVS6rQYY5i6aCcdmtTj4u7NrI6jPCg4MIDLzmjBV+MGMbD3ckqCEgm1d6bcFHt8Rl5nWhzdjDG5wJXAj0Bb4BaPplJKnZZfUrLYeiCPcUPaExCgrY3aIDMvk693TSYr+N8cCH0cm91GQnKCR1sdzhSOYBEJpqJwzDHGlAL+1zGilI8zxvDOop20alSXK3q1sDqOqiF/nhur4mu5Kfdoq8OZwvE+sAcIB5aKSBsg12OJlFKnZcWubJLTjjF2cHuCvGiwmPKsqubGWpG+wmPnPOXIcWPMW8BblTbtFZGh7ji5iAwH3gQCgWnGmEnHvR4KfAr0A7KB64wxe9xxbqX8zTsLd9IkIpSr+8ZYHUXVoKS7k2r8nKcsHCISCdwKxB63//3VObGIBAJTgAuBdGCtiMwxxmyptNsdwFFjTAcRuR74F3Bddc6rlD9at/coK1OzefrSrj4zw6ryXc7MVfUDsArYCNjdeO7+wE5jTCqAiHwBjAQqF46RwPOO778E3hERMf44+ESpapiyaOcfs6sq5WnOFI46xpiHPXDulkBapefpwFlV7WOMKRORHCAKOHz8wURkDDAGoHVr/cejao/N+3NYuDWLRy7sRFiIzluqPM+ZHrT/iMhdItJcRBr9/vB4MhcZYz4wxsQZY+Kio6OtjqNUjZm6eBcRoUHcOjDW6iiqlnDmzxMb8Aowgf+/DdcA7ap57gygVaXnMY5tJ9onXUSCgAZUdJIrpYBdh/L5YWMm9wxuT4O6wVbHUbWEM4XjEaCDMeYvl4eqaS3QUUTaUlEgrgduPG6fOcBtwErgb8BC7d9Q6v+9u3gXoUEB/P2ctlZHUbWIM4VjJ1Do7hM7+izGA/OouB33Y2PMZhGZCCQaY+YAH1FxqWwncISK4qKUAtKPFvJNUgY3n92GxvVCrY6jahFnCkcBkCwii4CS3zcaY6p1O67jGD9QcddW5W3PVvq+GLimuudRyh+9vyQVEbh7cHWvGivlGmcKxzeOh1LKS2TlFjMrMY2r+8bQvEFdq+OoWsaZkeOf1EQQpZTzPlq+m7JyO2MHt7c6iqqFnBk5vpsTTGpojNH2sVIWOFZo47NVe7m8VwtiG4dbHUfVQs5cqqq8IlQdKvocvG4ch1K1RcKveyiwlTNuSAero6ha6pQDAI0x2ZUeGcaYycClNZBNKXWc/JIypq/Yw0XdmtK5WYTVcVQt5cylqr6VngZQ0QLReQ2UssBnq/aSU1TKvUO1taGs40wBeK3S92VUrM1xrUfSKKWqVFxazrRluzm3Y2N6tYq0Oo6qxZy5q8ota28opapndmIah/NLuHdoH6ujqFrulH0cIvKAiNSXCtNEZL2IXFQT4ZRSFUrL7by/JJW4Ng05q63em6Ks5czsuH83xuQCF1ExpfktwKSTv0Up5U5fJ2WQcayIe8/vgIhYHUfVcs4Ujt9/Si8BPjXGbK60TSnlYek5+3nq28V0bhbGkE66ZICynjOFY52IzKeicMwTkQjcuxKgUuokxn/9EWWljQhvuFxbG8orOFM47gCeAM40xhQCIcDtHk2llAIqWhtrt0dRKmn8lDaJA/kHrI6klFMDAO3GmPXGmGOO59nGmN88H00pNWbWZwTb23As+DPKKSN+SbzVkZRyqsWhlLLAhoy9bEptS1FAIoUBv2Irt5GQnKCtDmU5LRxKeam7P/8ZIZAjwe/+cTtKuSnXVoeynDPjOF4Tke41EUYpVWHh1oMcyG7OsaAvKAs4+Md2W7mNFekrLEymlHNTjqQAH4hIEJAAzDTG5Hg2llK1V6GtjGe+2UzHJvWYe/8XhATNtjqSUn/iTOf4NGPMIOBWIBb4TUQ+FxGdikQpD3jzlx1kHCvipVE9CQnSq8nK+zj1UykigUAXx+MwsAF4WES+8GA2pWqdrQdy+WjZbq6La8WZsTq1iPJOzkyr/gZwOfAL8JIxZo3jpX+JyDZPhlOqNrHbDU99tZH6dYN5YkQXq+MoVSVn+jh+A542xhSc4LX+bs6jVK31xdo01u87xmvX9KJheIjVcZSqkjOFYwPQ+bipDnKAvdpJrpR7HMorYdKPKQxoF8Wovi2tjqPUSTlTOKYCfaloeQjQA9gMNBCRe4wx8z2YT6la4cW5WygutfPCVT10Pirl9ZzpHN8P9DHGxBlj+gF9gFTgQuDfp3NSEWkkIj+LyA7H14ZV7FcuIsmOx5zTOZdS3m75jsN8k7yfsUPa0z66ntVxlDolZwpHJ8dU6gAYY7YAXYwxqdU47xPAL8aYjlR0uj9RxX5FxpjejscV1TifUzLzMhk8fbBO6aBqTHFpOc98u4m2jcMZN6S91XGUcoozhWOLiLwrIoMdj6mObaFA6WmedyTwieP7T4ArT/M4bvXEvNf4de86ndJB1Zipi3ex+3AB8SN7UCc40Oo4SjnFmcJxG7ATeNDxSAVGU1E0TncQYFNjTKbj+wNA0yr2qyMiiSKySkQ8Wly2ZqWxKLE/kbY7dSI5VSN2HcrnvcW7uLJ3C87p2NjqOEo57aSd446Bf9OMMTcBr51gl/yTvHcB0OwEL02o/MQYY0TEVHGYNsaYDBFpBywUkY3GmF1VnG8MMAagdevWVcWq0ltrX6YguICI0msptW0gfkk8Uy6d4vJxlHKGMYYJX2+kTnAAEy7tZnUcpVxy0haHMaYcaCMiLt9UbowZZozpcYLHt8BBEWkO4PiaVcUxMhxfU4HFVHTMV3W+Dxwd+HHR0a4tr5mZl0lCcgJHAmdQIlupXzKWT9bP0VaH8piv1mewKvUIT4zoSnREqNVxlHKJM5eqUoFfReQZEXn490c1zzuHiktgOL5+e/wOItLQ0Y+CiDQGBgFbqnneE4pfGo/d2EHKORzyKkIAEUXjmbhY+zqU+x0tsPHiDyn0bR3J9We2sjqOUi5zpnDsAr537BtR6VEdk4ALRWQHMMzxHBGJE5Fpjn26AokisgFYBExy3NHldivTV2IrtwFQFnCAI8HvE2rvwaLN2lmp3G/Sj1vJKSrlxat6EhCgYzaU7znlAEBjzD8BRCTMseZ4tRljsoELTrA9EbjT8f0KoKc7zncqSXcnHZ+D8Z8nMW/zcDam59AzpkFNxFC1wJrdR5iVmMbd57Wja/P6VsdR6rQ4s5DTABHZAmx1PO/luCXXb4kIL17Vg+iIUB74IolCW5nVkZQfsJXZmfD1RlpG1uWBYR2tjqPUaXPmUtVk4GIgG8AYswE4z5OhvEFkWAivXduL3dkFxH+fYnUc5Qc+XJbKjqx8Jo7sTliIM7P9KOWdnFqPwxiTdtymcg9k8ToD2zfm7vPaM3PNPuZv1jus1Onbl13IW7/sYHj3ZlzQtaphS0r5BmcKR5qIDASMiASLyKNULCdbKzx8YSd6tKzPP/73G1m5xVbHUT7IGMMz324iKEB47gods6F8nzOFYyxwL9ASyAB6O57XCiFBAUy+rg9FpeU88t8N2O1VjVVU6sTmbsxkyfZDPHJRZ5o3qGt1HKWqzZk1xw8bY24yxjQ1xjQxxtzsuCuq1ujQpB7PXNaNZTsOk7Bij9VxlA/JLS7ln99toUfL+tw2MNbqOEq5hTNLx0YDdwGxlfc3xvzdc7G8z439W7No6yH+9eNWBraP0lsplVNenbeN7PwSProtjkAds6H8hDOXqr4FGgALgLmVHrWKiPCvq3vSICyYB75Iori0VtwfoKohOe0Y/1m1l1sHxHJGTKTVcZRyG2cKR5gx5h/GmNnGmP/9/vB4Mi8UVS+UV6/pxfaD+Uz6cavVcZQXKyu389RXG2kSEcojF3WyOo5SbuVM4fheRC7xeBIfMbhTNH8f1JbpK/awaOsJ52ZUiukr9rAlM5fnLu9ORJ1gq+Mo5VbOFI4HqCgexSKSKyJ5IpLr6WDe7PHhnenSLILHvtzA4fwSq+MoL7P/WBGv/7ydoZ2jGdHjRCsLKOXbnLmrKsIYE2CMqWOMqe94Xqt7husEB/Lm9X3ILS7j8S9/wxi9RVf9v+fnbMZuDBNH9kBEO8SV/3FmrioRkZtF5BnH81Yi0t/z0bxb52YRPDmiCwu3ZvHZ6n1Wx1Fe4uctB5m/5SAPXNCJVo3CrI6jlEc4c6lqKjAAuNHxPB/QpfGA0QNjGdwpmhe+38LOrDyr4yiLFZSU8dy3m+jcNII7z21rdRylPMaZwnGWMeZeoBjAGHMUcHlFQH8kIrxyzRmEhwZx/8xkSsr0Ft3abPKC7ezPKealUT0IDnRqGjilfJIzP92ljrXHDfwxINDu0VQ+pElEHf599RlsyczltfnbrY6jLLJlfy4f/7qHG/q3ol+bRlbHUcqjnCkcbwFfA01E5EVgOfCSR1P5mGHdmnLz2a35YGkqy3cctjqOqmHldsNTX28ksm4w/xjexeo4SnmcM3dVzQAeB14GMoErjTH/9XQwXzPhkm60jw7nkf8mc7TAZnUcVYM+X7OP5LRjPH1ZVyLD9Cqu8n/Orsex1RgzxRjzjjGm1kyp7oq6IRW36B4psPHEV3qLbm2RlVfMv3/ayqAOUVzZu6XVcZSqEdqD50Y9WjbgsYs7M2/zQWYnHr/2lfJH8d+nUFJqJ17HbKhaRAuHm915TjsGto/i+TlbSD2Ub3Uc5UFLtx/iuw37GTe0Pe2i61kdR6kao4XDzQIChNeu7UVIUAAPzUqmtFxvQPM3mXmZnPvxBTz5VTLtGodzz5D2VkdSqkZp4fCA5g3qMmlUTzak5/Dmgh1Wx1FuFr80nk2pLck4ZuOFq3oQGhRodSSlapQWDg8Z0bM518bFMGXxTlan1qoFE/1aZl4mn66bT/2yURQFLaFd0zKrIylV4ywpHCJyjYhsFhG7iMSdZL/hIrJNRHaKyBM1mdEdnru8O20ahfHw7A3kFJVaHUe5wcQl8UQUj8FOMbkh04lfEm91JKVqnFUtjk3AKGBpVTs4RqtPAUYA3YAbRKRbzcRzj/DQICZf34cDucU8/c0mvUXXx2XmZTJrbRqh9u4cDf6YYnOIhOQEDuQfsDqaUjXKksJhjEkxxmw7xW79gZ3GmFRjjA34Ahjp+XTu1btVJA8N68h3G/bzTXKG1XFUNUxY8G/qldxCccBvFAT+DEC5KddWh6p1vLmPoyVQeTBEumObz7lnSAfOjG3IM99sJu1IodVx1GlatrEZQjDZwe+AY8iGrdzGivQV1gZTqoYFeerAIrIAONHyZxOMMd964HxjgDEArVu3dvfhqyUwQHjjut6MmLyMB2clM2vM2QTp7Kk+5ectBykt7MHjF3fm3qHaclS1m8d+exljhhljepzg4WzRyABaVXoe49hW1fk+MMbEGWPioqOjqxPdI2IahvHCVT1Yt/coUxbtsjqOckFecSnPfFOxzsZd57azOo5SlvPmP3vXAh1FpK2IhADXA3MszlQtI3u35Ko+LXlr4Q7W7T1qdRzlpFfnbeNgXjEvX92TkCBv/iejVM2w6nbcq0QknYqVBeeKyDzH9hYi8gOAMaYMGA/MA1KA2caYzVbkdad/juxOs/p1uH9mIud+dKHekePl1u87yqer9nLr2W3o27qh1XGU8grij7eIxsXFmcTERKtjVGntniNc894KCgJ/4dqBNqZcqivxeiOXu7aIAAASIElEQVRbmZ3L315ObnEpPz88mHqhHusSVMpyIrLOGFPluLrKtN1tgZioEvJDviS8fBgzEpO01eGlPlyWyraDeUwc2UOLhlKVaOGwQPzSePKC/0up7Ce85Bb+ufgFqyOp46QeyufNX3ZwSc9mXNitqdVxlPIqWjhqWGZeJgnJCdjshRwNTiDY3ppZa9O01eFFjKlYCjY0KIDnL+9udRylvI4WjhoWvzQeu6mYar0oYCXFARsJt13Hs7+8bHEy9bv/JqazKvUIT47oSpP6dayOo5TX0cJRw1amr8RW7liTXOBo8EcEmgYs3RJqbTAFwKG8El78IYX+sY24/sxWp36DUrWQ9vjVsKS7k/6y7eFZyXy/cShpRwpp1SjMglTqdxO/30KRrZyXRvUgIECXglXqRLTF4QUevbgzArwy71TzPipPWrQ164+lYDs0ibA6jlJeSwuHF2gRWZcx57Vjzob9JO3TEeVWKCgp4+lvNtGhST1dClapU9DC4SXuHtyexvVCeWFuiq7bYYHXf95OxrEiXh7VU5eCVeoUtHB4iXqhQTx6USfW7T3KDxv11tyatCHtGAm/7uams1pzZmwjq+Mo5fW0cHiRa+Ja0aVZBJN+SqGkrNzqOLVCabmdJ77aSON6ofxjRBer4yjlE7RweJHAAGHCpV1JO1LEJyv2WB2nVvho+W5SMnOZOLI79esEWx1HKZ+ghcPLnNsxmiGdo3l74U6OFNisjuPX9mYXMHnBdi7s1pSLu59ozTGl1Ilo4fBCEy7pSqGtnDcXbLc6it8yxjDh600EBQQwcWR3RHTMhlLO0sLhhTo2jeD6M1vx2ep97MzKtzqOX/o6KYPlOw/z+PDONG9Q1+o4SvkULRxe6qELO1E3OJBJP6ZYHcXvZOeXEP/9Fvq2juTms9pYHUcpn6OFw0s1rhfKvUM7sCAlixU7D1sdx2WZeZkMnj7YK2f9fXFuCvklZbw86gydVkSp06CFw4vdPiiWlpF1eWFuCuV23xoU+M8lL7B833Lil8RbHeVPlm4/xFdJGYwd3J7OzXRaEaVOhxYOL1YnOJB/jOjClsxcvlqfbnUcp63dt4e5v55FdPEk/rNusde0Oops5Uz4ZiPtGodz79AOVsdRymdp4fByl5/RnN6tInll3jYKbWVWxzml7PwSbvt4DWJCCLY3p2HhJK6d9hVHveDW4skLtpN2pIiXRvWkTrBOK6LU6dLC4eVEhGcu60pWXgkfLE21Os5JFdnKufXjFRQUB5EVMpGMOmPJC/yOvQdiGPzqQj5btdeyS26bMnKYtnw318W14ux2UZZkUMpfaOHwAf3aNOLSns15f0kqB3KKrY5zQmXldu6bmcTm/QUcrTOZksAUjBRwNORDDtd9hKCQLJ7+ZhNXvLOcdXuP1Hi2J7/aSMOwEJ66pGuNnlspf6SFw0f8Y3gXyu2GV+d735odxhiem7OZBSkHqRv1PXmy7E+vF7ITe8M3ePuGPhwpsHH1uyt5eHYyWXk1UwSnr9jDxowcnru8Gw3CdFoRpapLVwD0Ea2jwhg9KJYPl6UyemAsPVo2sDrSH6Yu3sWM1fu4+7x2PHnJe8B7Ve57fpcmTFm0k2nLdjN/80EeHNaR2wbGEhzomb9h0o4U8tr87ZzfpQmXndHcI+dQqraxpMUhIteIyGYRsYtI3En22yMiG0UkWUQSazKjN7p3aAci6wbzohet2fHV+nRembeNK3q14B/DTz27bHhoEI8P78K8h84jLrYhL8xNYcSby/jVA2NVjDE8/c0mRCD+yh46rYhSbmLVpapNwChgqRP7DjXG9DbGVFlgaosGdYN5cFgnVqZm80tKltVxWLbjEI9/+RsD2kXxyjWuDaZr2zichNFnMu3WOGxldm6atppxM9aRcazIbfnmbNjPku2HePSizrSM1GlFlHIXSwqHMSbFGON9F+t9wI1ntaZddDgv/ZBCabndshyb9+dwz2fr6dCkHu/f2u+0Vs0TEYZ1a8r8h87jkQs7sXBrFhe8tpi3f9lBcWn11iM5Vmhj4ndb6BXTgNsGxlbrWEqpP/P2znEDzBeRdSIyxuow3iA4MICnRnQl9XABn6/eZ0mG9KOFjE5YS0SdIBJuP7Pa61jUCQ7kvgs6suDhwZzfpQmv/bydi95Yyi8pB0/7mC/OTeFYUSkvjzqDQJ1WRCm38ljhEJEFIrLpBI+RLhzmHGNMX2AEcK+InHeS840RkUQRSTx06FC183uzC7o2YWD7KCYv2E5OUWmNnvtYoY3RCWspLi1n+u393TqzbEzDMKbe1I8Zd55FSFAAd3ySyO0Ja9h9uMCl46zYeZj/rkvnrnPb0a1FfbflU0pV8FjhMMYMM8b0OMHjWxeOkeH4mgV8DfQ/yb4fGGPijDFx0dHR1f8AXkykYqXAY0WlTFm0s8bOW1xazphP17Evu5APbonz2FxPgzo05scHzuXpS7uyds9RLn5jKf/+aatTI+eLS8t56uuNtIkK48FhHT2ST6nazmsvVYlIuIhE/P49cBEVneoK6N6iAX/rG8P0X/ewL7vQ4+ez2w0Pz05mzZ4jvHptLwa09+zo6+DAAO48tx0LHxnMZb2aM3XxLi54bQnfbdh/0jvK3l64gz3Zhbx4pU4ropSnWHU77lUikg4MAOaKyDzH9hYi8oNjt6bAchHZAKwB5hpjfrIir7d69OLOBAYI//ppq0fPY4whfu4Wfth4gAmXdOWKXi08er7KmtSvw+vX9ubLsQNoGBbCfTOTuOHDVWw7kPeXfVMyc3l/SSpX943hnI6NayyjUrWNeMt4AHeKi4sziYm1Y9jH5AXbmbxgB1+OHUBcbCOPnGPaslRemJvC7YNiefaybpaNhyi3G2au2cer87eRV1zGrQPa8OCwThSWHea6L28gIv95Mo+VsuDhwTQKD7Eko1K+SkTWOTvsQQuHjyu0lTH01cU0a1CXr+8Z6PaFib7bsJ/7ZiYxokcz3rmxr1fcoXS0wMar87fx+Zp9RIWH0Lz5WpbuWU3D0jFMvq43V/ZpaXVEpXyOK4XDa/s4lHPCQoJ49KLObEg7xne/7XfrsVfuyuaR2Rs4M7Yhb1zX2yuKBkDD8BBevKon340/h2YNgtm4oxcNS8dQEpjMWR30R1opT9N/ZX7g6r4xdG9Rn3//tK3aA+d+t+1AHmP+k0jrqDA+vDXOKzuae7RsQKu2szkW+hZFAevJCX2XF5a+YHUspfyeFg4/EBBQcXtuxrEiPv51d7WPdyCnmNEJa6gTHMj0288kMsw7+wsy8zKZviGBnID5ZIU+S5HJICE5wWtWHFTKX2nh8BMD2zdmWNemTF20i8P5Jad9nNziUkYnrCG3qJTpt59JTMMwN6Z0r/il8djNn6ddKTflXrfOuVL+RguHH3nyki4Ul5bzxs/bT+v9tjI7Y/+zjp1Z+bx3Sz+6t/CeqdtPZGX6Smzlf16S1lZuY0X6CosSKVU76HocfqR9dD1uPrsNn67cw20DY+nU1PmR3Xa74bEvN7BiVzavXdOLczt6/+j7pLuTrI6gVK2kLQ4/c/8FHQkPDeKlH1Jcet+/5m3l2+T9PHZxZ67uF+OhdEopf6CFw880Cg/h/vM7snjbIZZud26yx09W7OH9JancdFZrxg1p7+GESilfp4XDD906sA2tG4Xx4twUyu0nH+D506YDPP/dZoZ1bcI/r+iuq+QppU5JC4cfCg0K5IkRXdh2MI/ZiWlV7rdu7xEe+CKJXjGRvH1DX4I8tO63Usq/6G8KPzWiRzPi2jTktfnbyS/563Tkuw7lc8cniTRvUIePboujboj3DfBTSnknLRx+6vc1Ow7nl/De4l1/ei0rr5jbPl5DoAif/L0/UfVCLUqplPJFWjj8WJ/WDRnZuwUfLktl/7EiAPJLyvj79LVk59v4ePSZtIkKtzilUsrX6DgOP/fYxZ35cdMBJn6fzJbSp2hROpGUzDw+vLUfvVpFWh1PKeWDtHD4uZiGYdx5TlumLt5FUcD57LPnMmlUT87v0tTqaEopH6WXqmqBUXH1KJdj1LX3pSDkvwzp5p2TFiqlfIMWjlrg9dUvcSz0DY4GfURu8EydBFApVS1aOPxcZl4mCckJ5Ms6coO/xma36dTjSqlq0cLh53TqcaWUu2nh8HM69bhSyt30rio/p1OPK6XcTVscSimlXKKFQymllEu0cCillHKJFg6llFIu0cKhlFLKJWLMyVeI80UicgjY64FDNwYOe+C4VtLP5Bv0M3k/X/88bYwx0c7s6JeFw1NEJNEYE2d1DnfSz+Qb9DN5P3/7PCejl6qUUkq5RAuHUkopl2jhcM0HVgfwAP1MvkE/k/fzt89TJe3jUEop5RJtcSillHKJFg4XiUi8iPwmIskiMl9EWlidqbpE5BUR2er4XF+LiM8vRi4i14jIZhGxi4jP3ukiIsNFZJuI7BSRJ6zO4w4i8rGIZInIJquzuIOItBKRRSKyxfEz94DVmTxNC4frXjHGnGGM6Q18DzxrdSA3+BnoYYw5A9gOPGlxHnfYBIwCllod5HSJSCAwBRgBdANuEJFu1qZyi+nAcKtDuFEZ8IgxphtwNnCvn/x/qpIWDhcZY3IrPQ0HfL6TyBgz3xhT5ni6CoixMo87GGNSjDHbrM5RTf2BncaYVGOMDfgCGGlxpmozxiwFjlidw12MMZnGmPWO7/OAFKCltak8S9fjOA0i8iJwK5ADDLU4jrv9HZhldQgFVPzySav0PB04y6IsygkiEgv0AVZbm8SztHCcgIgsAJqd4KUJxphvjTETgAki8iQwHniuRgOehlN9Jsc+E6hods+oyWyny5nPpFRNEZF6wP+AB4+7MuF3tHCcgDFmmJO7zgB+wAcKx6k+k4iMBi4DLjA+co+2C/+ffFUG0KrS8xjHNuVlRCSYiqIxwxjzldV5PE37OFwkIh0rPR0JbLUqi7uIyHDgceAKY0yh1XnUH9YCHUWkrYiEANcDcyzOpI4jIgJ8BKQYY163Ok9N0AGALhKR/wGdATsVM/CONcb49F+BIrITCAWyHZtWGWPGWhip2kTkKuBtIBo4BiQbYy62NpXrROQSYDIQCHxsjHnR4kjVJiIzgSFUzCZ7EHjOGPORpaGqQUTOAZYBG6n4vQDwlDHmB+tSeZYWDqWUUi7RS1VKKaVcooVDKaWUS7RwKKWUcokWDqWUUi7RwqGUUsolWjhUrSMikSIyrtLzISLyvYvHGO3pmZFF5HkRedST51DqdGjhULVRJDDulHud3GjA56fUV+p0aOFQtdEkoL1jTZVXHNvqiciXjnVJZjhGAyMi/URkiYisE5F5ItJcRP4GxAEzHMeoKyLPishaEdkkIh/8/v7fiUgDEdkrIgGO5+EikiYiwSJyl+O9G0TkfyISdnxgEVn8+7oiItJYRPY4vg90rKey1rGeyt2O7c1FZKkj3yYROddD/y1VLaSFQ9VGTwC7jDG9jTGPObb1AR6kYt2LdsAgx/xDbwN/M8b0Az4GXjTGfAkkAjc5jlEEvGOMOdMY0wOoS8W8X38wxuQAycBgx6bLgHnGmFLgK8d7e1ExJfcdLnyWO4AcY8yZwJnAXSLSFrjRcfzeQC/HuZVyC53kUKkKa4wx6QAikgzEUjFVSQ/gZ0cDIhDIrOL9Q0XkcSAMaARsBr47bp9ZwHXAIirmnZrq2N5DRF6g4hJaPWCeC7kvAs5wtIIAGgAdqZjn6mNH8fvGGKOFQ7mNFg6lKpRU+r6cin8bAmw2xgw42RtFpA4VRSDOGJMmIs8DdU6w6xzgJRFpBPQDFjq2TweuNMZscMxSPOQE7y3j/68QVD62APcZY/5SbETkPOBSYLqIvG6M+fRkn0MpZ+mlKlUb5QERTuy3DYgWkQFQMXW2iHQ/wTF+/0V+2LEmw984AWNMPhUtgTeB740x5Y6XIoBMR+vgpiqy7KGi2HDc8ecB9zjei4h0cvSftAEOGmM+BKYBfZ34vEo5RVscqtYxxmSLyK8isgn4EZhbxX42xyWgt0SkARX/XiZTcRlqOvCeiBQBA4APqVjn/AAVxaEqs4D/8udWxTNUrBh3yPH1REXtVWC2iIw5Lu80Ki6rrXd0yB8CrnQc/zERKQXyqVixUim30NlxlVJKuUQvVSmllHKJFg6llFIu0cKhlFLKJVo4lFJKuUQLh1JKKZdo4VBKKeUSLRxKKaVcooVDKaWUS/4PNtJ4RNsWkIQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(thetas, diff_phys_sim, 'g^')\n",
    "plt.plot(thetas, diff_phys_sim)\n",
    "plt.xlabel('theta values')\n",
    "plt.ylabel('energy sums')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
